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ABSTRACT 


Most  self-propelled  vessels  moving  on,  or  under,  the  ocean  surface, 
contain  rotating  machinery  that  radiate  finite  bandwidth  signals  into  the  water. 
Empirical  evidence  suggests  that  the  signal  bandwidth  estimated  with  a  far  field 
receiver  is  often  greater  than  expected.  This  thesis  investigates  the  use  of  an 
acoustic  propagation  model  to  predict  the  received  bandwidth  of  sinusoidal 
signals  when  both  the  source  and  the  receiver  are  in  motion.  The  bandwidth 
parameter  is  calculated  from  the  multi-frequency  transmission  loss  (TL)  predicted 
with  a  re-written  version  of  K.  Smith’s  Monterey-Miami  Parabolic  Equation 
(MMPE)  model,  including  both  receiver  and  source  motion.  The  results  for 
various  propagation  environments  allow  exploration  of  the  characteristics  of 
received  bandwidth,  predicted  from  sources  on  the  surface  or  at  depth.  The 
dependency  of  aggregate  bandwidth  upon  conditional  parameters  such  as  range, 
depth,  and  normalized  pressure  are  also  evaluated.  In  addition  to  modeling 
results,  this  thesis  documents  a  new  implementation  of  the  MMPE  model  for 
narrowband  signals  using  only  the  MATLAB  programming  language.  A  MATLAB 
version  has  the  inherent  advantages  of  increased  flexibility  and  portability.  A 
MATLAB  implementation  of  a  range  dependent  ray  trace  function  based  upon  a 
Runge-Kutta  integration  of  the  eikonal  equations  is  also  presented. 
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I.  INTRODUCTION 


This  thesis  is  motivated  by  several  remarks  in  the  literature  concerning 
spreading  of  a  continuous  wave  (CW)  signal  in  the  frequency  domain  after 
propagation  in  an  ocean  waveguide  (see  refs.  1,  2,  and  3  for  example).  In  many 
cases,  the  frequency  spreading  can  be  attributed  to  differences  in  the  received 
Doppler  (i.e.  differential  Doppler)  of  various  multipath  arrivals  when  the  source 
and/or  receiver  are  in  motion.  Experiments  also  suggest  that  a  Doppler  shift 
imparted  on  transmitted  signals  due  to  surface  motion  is  another  factor  (see  refs. 
4  and  5  for  example)  that  can  cause  frequency  spreading  even  when  the  source 
and  receiver  are  stationary.  The  spread  of  power  across  frequency  can  be 
quantified  into  a  single  bandwidth  parameter  that  is  useful  for  aggregate  signal 
analysis  or  source  classification.  In  the  case  of  source/receiver  motion,  it  is  likely 
that  signal  bandwidth  is  dependent  upon  the  parameters  defining  the  acoustic 
propagation  channel  and  the  relative  positions  and  motion  of  the  source  and 
receiver.  To  predict  the  dependencies  due  to  differential  Doppler,  a  model  of  the 
broadband  transmission  loss  (TL)  in  the  presence  of  source  and  receiver  motion 
is  required. 

The  literature  related  to  acoustic  propagation  models  incorporating  the  effects  of 
source  and  receiver  motion  is  diverse.  Much  of  the  published  work  focuses  upon 
approaches  based  upon  ray  theory.  A  series  of  papers  by  Jacobson  and 
colleagues®’^’®’®’^°”^^  use  ray  theory  to  evaluate  the  effect  of  source  motion  on 
transmission  loss  for  specific  sound  speed  profiles  (e.g.,  bilinear  or  iso-speed 
SSP).  A  total,  time-dependent,  pressure  field  amplitude  and  phase  is  computed 
by  summing  individual  ray  contributions,  where  parameters  such  as  travel  time, 
received  angle,  and  frequency  are  modified  by  relative  motion.  The  papers  differ 
in  their  focus  on  specific  effects  associated  with  short,  intermediate  and  long 
ranges,  as  well  as  various  depths  and  SSPs.  A  pair  of  papers  by  Flanagan  and 
Weinberg^^’^^  also  use  ray  theory  to  look  at  the  coherence  length  and  bandwidth 
of  multiple  CW  signals  in  the  presence  of  radial  and  skewed  source  motions. 
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Periodic  source  motion  is  shown  to  cause  an  amplitude  modulation  that 
“transfers  power”  into  sidebands,  thereby  increasing  the  signal’s  effective 
spectral  bandwidth. 

References  14  &  15  are  among  the  first  papers  to  incorporate  source 
and/or  receiver  motion  into  a  normal  mode  model.  Hawker^'^  derived  an 
approximate  expression  for  the  time-dependent  pressure  field  as  a  function  of 
wavenumbers  (eigenvalues),  relative  source  speed,  mode  group  velocities,  and 
mode  (eigen)  functions.  Included  in  the  derivation  is  an  approximate  formulation 
for  the  Doppler  shifted  frequencies  of  the  individual  modes  as  a  direct  function  of 
the  phase  and  group  velocities.  Another  formulation  of  the  pressure  field 
associated  with  a  moving  source  was  independently  derived  in  an  earlier 
publication  by  Neubert^^.  The  unique  approximations  and  assumptions  invoked  in 
this  work  resulted  in  a  somewhat  different  formulation  than  that  presented  by 
Hawker.  Neubert,  however,  appears  to  focus  on  how  a  conventional  stationary 
source/receiver  model  can  be  modified  to  estimate  the  field  for  moving  source  or 
receivers. 

Hawker’s  approach  was  referenced  many  years  later  by  Schmidt  and 
Kuperman,^^  where  it  was  compared  to  another  model  derived  from  a  wave- 
number  integration  approach.  The  later  formulation  was  used  to  predict  the 
received  frequency  response  of  a  band-pass  signal  emitted  from  a  moving 
source  in  a  shallow  water  waveguide.  Song  and  Baggeroer^^  also  referred  to 
Hawker’s  results  in  their  development  of  an  algorithm  to  estimate  source  velocity 
by  computing  the  Doppler  shifts  of  individual  modes. 

Methods  for  including  source  or  receiver  motion  in  a  parabolic  equation 
(PE)  based  acoustic  propagation  model  have  also  been  explored.  Howell, 
Jacobson  and  Seigman^®  published  a  solution  that  applied  a  Galilean 
transformation  to  the  wave  equation.  This  equation  is  solved  using  a  time 
harmonic  pressure  solution  containing  a  frequency  term  that  is  modified  by 
source  and  receiver  motion.  The  Helmholtz  equation  is  converted  to  a  parabolic 
equation  using  the  standard  far-field  approximations.  In  turn,  this  is  converted  to 
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the  narrow  angle  parabolic  equation  by  approximating  a  complicated  expression 
involving  the  refraction  index  by  a  simpler  term  containing  an  “effective”  refraction 
index.  The  result  is  then  solved  with  standard  marching  algorithms.  Numerical 
results  were  presented  for  a  waveguide  with  a  constant  speed  of  sound,  although 
the  technique  appears  sufficiently  general  to  be  applied  to  any  spatially  variant 
sound  speed  profile. 

Another  approach  is  considered  by  Smith^^  who  represents  source  motion 
by  modifying  the  starting  field  that  initializes  a  Split-Step  Fourier  (SSF)  PE 
marching  algorithm. Since  the  source  frequency  of  the  starting  field  is  a 
function  of  the  transmitting  angle  and  source  speed  (relative  to  the  environment), 
the  transmitting  wavenumbers  also  experience  a  similar  adjustment.  Thus,  a  new 
starting  field  can  be  calculated  that  incorporates  the  frequency  dependent  source 
distribution  as  a  function  of  vertical  wavenumbers.  Smith  also  describes  a 
method  for  representing  receiver  motion  by  numerically  interpolating  the  received 
pressure  field  across  the  frequency  grid  at  each  wavenumber  and  range.  At  a 
particular  range,  multiple  responses  in  the  vertical  wavenumber  domain  versus 
frequency  are  “re-mapped”  to  a  new  wavenumber-frequency  domain.  This 
method  lends  itself  to  a  very  straightforward  implementation. 

A  fair  amount  has  also  been  written  regarding  the  effects  of  a  time  varying 
ocean  surface  on  the  characteristics  of  a  CW  signal. Although,  these 
effects  are  sometimes  important,  this  thesis  will  concentrate  on  the  effects  of 
source/receiver  motion  when  the  surface  is  flat. 

Given  the  options  identified  above,  we  have  chosen  to  base  our  study  on 
the  PE  acoustic  propagation  model  implemented  by  Smith.  To  facilitate  upgrades 
and  numerical  experimentation,  the  model  was  re-written  in  the  MATLAB 
programming  language  (vice  the  original  FORTRAN  code).  The  implementation 
includes  the  Doppler  effects  on  a  sinusoidal  signal  from  a  moving  source,  and  the 
frequency  shift  caused  by  a  moving  receiver.  The  resulting  code  has  been 
named  the  MATLAB  Monterey-Miami  Parabolic  Equation  (M3PE)  model.  The 
theoretical  background  of  a  PE  TL  model,  including  the  approach  for 
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implementing  source  and  receiver  motion  is  described  in  Sections  II  and  III.  In 
Section  IV,  the  functionality  of  the  new  implementation  is  validated  by  comparing 
transmission  loss  results  against  outputs  also  generated  by  the  original  MMPE 
and  a  code  based  upon  a  finite  difference  approach  for  several  benchmark 
environments. 

The  M3PE  model  was  used  to  conduct  a  series  of  numerical  experiments 
intended  to  explore  the  effective  bandwidth  of  the  TL  for  various  types  of  ocean 
environments  and  source  depths.  This  bandwidth  parameter  is  defined  in  Section 
V.  Numerical  results  provided  in  Section  VI  indicate  the  possibility  for  several 
trends.  First  it  was  apparent  that  the  calculated  bandwidth  parameter  often 
follows  the  level  of  TL  in  range  and  depth  space.  That  is,  higher  TL  often  results 
in  higher  bandwidth.  Second,  the  bandwidth  of  TL  is  influenced  by  the  tendency 
of  the  waveguide  to  propagate  high  angles.  Higher  bandwidths  were  observed 
when  acoustic  conditions  allowed  the  propagation  of  significant  energy  at  high 
angles,  as  opposed  to  smaller  bandwidths  that  resulted  from  conditions  that 
stripped  out  high  angle  modes.  Finally,  for  the  ocean  conditions  evaluated,  there 
was  a  tendency  for  sources  near  the  surface  to  result  in  wider  TL  bandwidths 
than  sources  positioned  deeper  in  the  water  column  and  closer  to  the  sound  axis. 
These  trends  require  further  study  since  it  is  not  completely  apparent  how 
observable  they  are  in  real  oceans. 

A  user’s  guide  for  operation  of  the  M3PE  code  is  provided  in  Appendix  A. 
In  addition  to  transmission  loss  and  bandwidth  assessment,  the  code  also 
includes  a  MATLAB  implementation  of  a  range-dependent  ray  trace  function.  The 
theoretical  background  and  implementation  details  associated  with  this  function 
are  described  in  Appendix  B. 
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II.  THEORETICAL  BACKGROUND 


A.  DERIVATION  OF  THE  SPLIT-STEP  FOURIER  PE  ALGORITHM 

The  following  description  closely  follows  discussions  presented  by 
Smith^°,  and  Thomson  and  Chapman^^  Much  of  Smith’s  work  is  based  upon  the 
approach  originally  published  by  Tappert^^’^^. 

If  T’(r,z,®t)  =  /»(r,z)e is  the  time-harmonic  solution  to  the  two- 
dimensional,  angularly  symmetric  wave  equation,  then  its  homogenous  solution 
in  cylindrical  coordinates  satisfies  a  Helmholtz  equation  of  the  form. 


iA 

r  dr 


dp^  d^  p 


dr 


+  +  {r,z)p  =  0  , 

dz 


(1) 


where  co  is  the  angular  frequency,  ko  =  (olco  is  the  reference  wave  number,  co  is 
the  reference  speed  of  sound,  and  n(r,z)=  co/c(r,z)  is  the  range  and  depth 
dependent  index  of  refraction.  This  describes  the  pressure  at  r  >  0  due  to  a  point 
source  at  r=0  and  some  variable  depth  for  a  constant  density  waveguide. 

The  cylindrical  spreading  term  can  be  eliminated  (for  convenience)  by 
substituting  p{r,z)  =  ^u{r,z) ,  yielding 


d^U  d^U  ,22/  X  A 

^  +  ^  +  {r,z)u  =  0  . 
dr  dz 

If  we  substitute  the  following  operators. 


d 

^op=^  and 


Qop  -Jn^  + 


1 

kl  dz^ 


Eq.  (2)  can  be  rewritten  as 

This  can  subsequently  be  factored  into 


(2) 


(3) 


(4) 
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[Pop  +  \Pop  -  iKQop  )m  =  0 


(5) 


where  the  commutator  operator  [Pop,Qop]  =  PopQop  -  QopPop  has  been  assumed  to 
equal  zero,  as  is  true  for  layered  media^^  Equation  (5)  represents  the 
factorization  of  the  field  into  incoming  and  outgoing  waves,  where  the  outgoing 
wave  must  satisfy 

PopH=ik,Q,pU.  (6) 

This  treatment  can  be  restricted  to  outgoing  waves  for  environments  where  the 
backscattered  field  is  small  relative  to  the  outgoing  field. 

The  field  function  u  for  the  outgoing  wave  can  now  be  decomposed  into 
two  parts  associated  with  a  slowly  varying  envelope  and  an  oscillating  phase 
function  as  defined  by 

(7) 

When  this  is  substituted  into  Eq.  (6),  the  result  is  the  following  expression  for  the 
outgoing  field, 

^  =  ,  (8) 

dr 


where  the  PE  field  term  y/,  is  related  to  acoustic  pressure  by  the  equation 


p{r,z) 


Vv 


W{r,z)e 


ik„r 


(9) 


Equation  (8)  is  a  differential  equation  of  the  parabolic  form,  vice  the 
elliptical  form  of  Eq.  (4).  To  solve,  though,  we  need  an  approximate  expression 
for  the  Qop  operator  that  facilitates  solution  of  a  differential  equation  with  a 
sequential  marching  algorithm.  In  our  case,  the  algorithm  propagates  the  solution 
in  range  using  the  representation 

y/(r  +  Ar)  =  (b(r)y/(r),  (10) 


where 
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Uop=\-^JY+s  =\-n  =  \--^  .  (17) 

c(r,z) 

Equation  (16)  is  called  the  kinetic  energy  (KE)  operator  and  Eq.  (17)  is 
called  the  potential  energy  (PE)  operator.  Using  Eqs.  (10)  and  (11),  we  can 
implement  the  marching  algorithm  using  the  substitution 

(//(r  +  Ar)  =  Q>{r)xi/{r)  =  _  (18) 

This  manipulation  has  effectively  split  the  original  operator  Qop,  into  an 
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operator  in  the  z  domain  {Uop),  and  an  operator  in  wavenumber  space  {Top). 
Next,  we  invoke  the  approximation,  by  using  the  Campbell- 

Baker-Hausdorff  expansion. Since  Top  and  Uop  are  small  the  commutator 
term  can  also  be  neglected.  (We  note  that  the  approximation  is  not  necessarily 
valid  for  operators  as  it  is  for  algebraic  variables,  but  it  is  assumed  to  be 
sufficiently  accurate  here.)  The  resulting  propagation  algorithm  becomes 


\l/{r  +  Ar)  =  e 


Ar 

-ikQ — t/op(r+Ar)  -j  x  t 

2 


-ik„—Uop(r) 

e  2  y/{r)  , 


(19) 


where  the  potential  energy  operator  was  segmented  to  center  the  result  at  each 
range  step. 

To  apply  the  KE  operator  in  practice,  we  require  a  scalar  version  of  Top 
defined  as^° 


ro,(*,)=i 


(20) 


fop{k,)  is  multiplied  by  the  transformed  field  in  wavenumber  space.  Appendix  C 
explains  how  this  is  equivalent  to  applying  the  original  operator . 


Then,  Eq.  (19)  is  implemented  by  the  following  sequence  of  Split-Step 
Fourier  (SSF)  transforms  and  propagation  functions^°: 


y/{r  +  Ar,z)  =  ^  ^  ^^ikM^p(K) 


xDFT 


Ar 

-ik„—Uop{r,z) 


xy/{r,z) 


,  (21) 


where  we  assume  the  MATLAB  convention  for  the  Discrete  Fourier  Transform, 
i.e. 

DFT-^  kv{k)  =  "^if,{z)e~‘'^'^  ,  k  =  0,l,...,N-l  ,  (22) 

z=0 

1  i—kz 

IDFT-^  y/{z)  =  -Y,^{ky  ^  ,  z  =  0,l,...,iV-l  .  (23) 

A  k=o 

Equation  (21)  concisely  represents  a  series  of  steps.  First  the  original  field 
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function  y/{r,z)  is  multiplied  by  the  PE  operator  (evaluated  at  range  r)  in  the 
depth  domain.  The  product  is  Fourier  transformed  into  wavenumber  space  using 
Eq.  (22).  Then  the  scalar  equivalent  of  the  KE  operator  is  applied  via 
multiplication  prior  to  inverse  transforming  the  product  back  into  the  depth 
domain.  Finally,  the  result  is  multiplied  by  the  PE  operator  evaluated  at  range 
r+Ar. 


In  practice,  the  variables  z  and  r  are  replaced  by  the  discrete  values  z„  and 
rj.  First  r,  is  defined  as 

rj=[j-^)^r  ,  7  =  1,2,3,...  (24) 


and  Zn  is  defined  on  the  half  space  mesh  as 


r 


n - 

V  2y 

1 


Az, 


N-n- 
V  2y 


N. 


iSz,  n=-^  +  \,...,N^ 


(25) 


For  conciseness,  we  will  drop  the  “z”  subscript  on  the  wavenumber, 
k^^k,  and  note  that  all  further  references  to  wavenumber  imply  the  vertical 
component,  unless  specifically  stated  otherwise.  The  sampling  in  wavenumber 
space,  kn,  is  then  defined  as 


k„  = 


ntsk ,  n  =  l,..., — ^ 

2 

-{N,-n)^k,  n  =  ^  +  \,...,N^ 


where 


M  ^  .  *2/) 

Ak  =  —  and  Az  =  — 

D  N 


(26) 


(27) 


and  D  is  the  maximum  depth  of  the  water  column. 

Up  to  this  point,  the  fluid  density  has  been  assumed  to  be  constant. 
However,  a  depth  variant  density  can  be  accommodated  in  the  solution  to  the 
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wave  equation  by  replacing  the  index  of  refraction,  n,  in  Eq.  (2),  with  an  effective 
index  given  by 


/  \  2 

ivV-^  -Vp  .  (28) 

2^0  P  ^\P  ) 

According  to  Smith, this  can  be  represented  in  the  marching  algorithm  by 
adding  a  second  term  to  the  potential  energy  propagator  as 

U{z)  =  U,{z)  +  U,{z),  (29) 

where  U\{z)  is  defined  by  Eq.  (17),  and  Uiiz)  is  defined  by  the  approximation 


(30) 


where  H”{z-z^)  is  the  second  derivative  of  a  density  transition  function  to  be 
defined  later  in  Section  lll-B. 


B.  THE  SOURCE  FUNCTION  INCLUDING  SOURCE  MOTION 

The  PE  marching  algorithm  requires  an  initial  condition  at  r  =  1  m  for  all  z. 
As  described  by  Smith,^°  an  omni  directional  point  source  can  be  specified  by  a 
pair  of  delta  functions  for  the  source  at  depth  z^,  and  its  image  above  the  surface. 
A  Fourier  transformation  converts  this  representation  back  to  the  wavenumber 
domain,  resulting  in  the  expression 

y/{r  =  Q,k^)  =  a\  [^(z-zJ-^(z  +  zJ]e^'''"^Jz  = -2zasin(A:„zJ,  (31) 

J-oo 

where  a  is  the  scale  factor  defined  as 


and  Ro  is  the  reference  distance  of  1  m. 

Next,  the  source  function  is  modified  by  a  function  that  corrects  far  field 
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errors  in  the  PE  solution.  This  function,  suggested  by  Thomson  and  Bohun,^°  is 
defined  as 


F{kJ 


1 

V 


-1/4 


(33) 


The  correction  function  is  set  to  zero  when  \kn\  >  ko,  since  this  represents  an 
evanescent  mode.  An  illustration  of  this  function  for  a  particular  value  of  ko  is 
shown  in  Fig.  1  below.  It  is  observed  that  the  correction  factor  in  the  starting  field 
is  most  significant  at  the  higher  wavenumbers  (i.e.,  larger  angles). 


WAv(»numb«r  for  sourco  fiirtdion 


Figure  1 .  Wavenumber  taper  for  source  function 

ik  — 

Finally,  an  extra  phase  term  of  e  "  ^  is  applied  to  account  for  the  half  cell 
depth  grid  defined  in  Eq.  (25).  The  result  is  a  wide  angle  point  source,  a.k.a. 
starting  function,  specified  in  the  wavenumber  domain  as 


^(r  =  0,^„)=-2z 


'Ijlkn 


^0  y 


-1/4 


Sin' 


ik^  — 
,  2 


(34) 


The  starting  function  presented  thus  far  is  only  valid  for  a  stationary  omni¬ 
directional  source.  To  account  for  source  motion,  we  recall  that  A:  =  A/g sin(6') , 
which  implies  that 


=sin' 


yko  j 


(35) 
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where  On  is  the  discretely  sampled  angle  from  the  horizontal  (positive  down) 
associated  with  each  wavenumber  within  the  bounds  \kn\  <  h.  Since  the  source 
function  is  defined  for  a  large  sector  of  wavenumber  space,  Eq.  (35)  implies  that 
the  source  is  also  generating  energy  at  many  vertical  angles. 

It  is  also  well  known  (see  [31],  for  example)  that  the  frequency  of  a 
continuous  wave  source  moving  along  the  horizontal  axis,  and  received  on  the 
same  horizontal  axis,  is  related  to  the  transmit  frequency  by  the  expression 


/'=/r 


(36) 


where  fr  is  the  original  transmit  frequency,  c  is  the  speed  of  sound,  is  the 
source  speed,  and  Vr  is  the  receiver  speed  (where  the  sign  of  and  Vr  infers  the 
horizontal  direction  of  motion).  If  the  receiver  is  stationary,  Eq.  (35)  can  be 
rewritten  as 


f'=fr 


1-^ 


^  V  ^ 
1  +  ^ 


V 


c  J 


(37) 


where  the  last  approximation  is  based  upon  a  binomial  series  expansion.  If  the 
source  is  also  moving  at  an  angle  (j)s  to  the  horizontal,  and  we  are  interested  in 
the  component  of  transmission  along  the  angle  0  relative  to  horizontal,  then  the 
general  expression  for  the  frequency  is 


f'-fr 


f  V  ^ 

1  +  — cos(6*-^4) 

VC  j 


(38) 


This  can  be  generalized  for  all  transmission  angles  On  by  the  expression 


fn'-fr  1  +  — COs(^„ 

V  c 


(39) 


Equations  (35)  and  (39)  define  a  relationship  between  the  transmitted 
frequency  and  the  vertical  direction  of  transmission.  Thus,  even  a  single 


12 


frequency  sinusoidal  source  will  produce  energy  across  a  bandwidth  of 
frequencies  as  a  function  of  direction  in  the  vertical  plane  when  the  source  is  in 
motion.  An  example  of  the  spread  of  frequencies  versus  wavenumbers  and 
associated  angles  for  a  100  Hz  source  traveling  at  20  knots  is  illustrated  in  Fig.  2 
for  three  different  source  traveling  directions,  (ps. 


Source  freguencv  vs  propegattng  wavenumbers 


Source  iretuency  vs  propapabng  ckrecnon 


(a)  (b) 

Figure  2.  Example  frequency  response  vs  wavenumber  (a)  or  angle  (b)  for  a  20 
knot  source  moving  at  the  grazing  angles  0°  (blue),  45°  (red),  or  -45° 

(green)  from  the  horizontal 


Numerical  implementation  of  this  effect  requires  distribution  of  the  source 
function  defined  in  (34)  onto  a  discrete  grid  having  coordinates  of  wavenumber 
versus  frequency.  We  performed  this  distribution  using  a  nearest  neighbor  rule 
followed  by  linear  interpolation.  The  process  begins  by  the  following  discrete 
assignments: 


j  a-y/{r  =  0,kj  ,  /;.</„'</, 

[(l  -  a)  ■  t//{r  =  0,  )  ,  otherwise 


(40) 


where  fi  represents  the  discrete  frequency  grid  for  which  the  PE  field  is 
evaluated,  and  f„’  is  defined  by  Eq.  (39).  Also  a  is  a  mixing  parameter  that  was 
nominally  set  to  0.9995  in  the  M3PE  implementation  to  provide  approximately  66 
dB  of  cell  attenuation. 

The  grid  frequencies  are  defined  as 
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(41) 


N.. 


/,=/r-^A/  +  (/-l)A/  ,  i  =  \,2,...,N. 


where 


A/ =  J/iVy  . 

(42) 

/max  =max(/„')  , 

(43) 

/min  =min(/„')  Vn  =  l...iV,  . 

(44) 

Here,  A^is  the  number  of  frequency  cells  and  n=l,2,...,Nz  defines  the  index  of 

Each  cell  is  then  smoothed  across  two  adjacent  frequency  bins  using 
linear  interpolation.  That  is,  for  each  bin  where/]  < /„'</]+i,  the  complex  source 

amplitudes  y/'[r  =  Q,k„,f)  and  =  and  y/'[r  =  Q,k^,f,_y)  are 

modified  as  follows: 


W 


(/:-,)= 


vif:) 

4/ 

.  v\f. 

V) 

.  f.>f: 
.  /i/.' 


vif.hvif.) 


1- 


(/„’-/■) 


A/ 


W 


fi^fn 

,  f,>fn 


(45) 


(46) 


(47) 


where  the  dependencies  upon  r  and  kn  are  dropped  for  simplicity. 

It  is  also  noted  that  any  constraints  on  the  directional  response,  as 
required  (for  example)  to  represent  transmission  from  a  phased  array,  can  be 
accommodated  by  simply  scaling  the  source  function  along  the  kn  axis  by  a 
weight  vector  that  is  the  inverse  Fourier  Transform  of  the  desired  vertical  beam 
pattern. 
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C.  INCORPORATING  RECEIVER  MOTION 

As  pointed  out  by  SmithJ®  the  influence  of  receiver  motion  cannot  be 
represented  by  a  simple  one-time  reshuffling  of  the  start-up  function  across  a 
frequency  grid,  as  described  for  the  source  motion  case  (see  Eqs.  44-46).  This  is 
because  the  effects  of  receiver  motion  will  vary  at  each  range  step,  and  because 
an  analytical  form  of  the  pressure  response  is  no  longer  available  after  the 
source  has  propagated  through  the  medium.  Any  realization  of  the  frequency- 
dependent  transmission  loss  at  a  particular  range  and  depth  must  wait  until  the 
source  function  propagates  through  the  medium.  Since  the  propagated  field  is 
only  defined  numerically,  receiver  motion  can  also  only  be  represented  with  a 
numerical  approach. 

The  goal  of  the  algorithm  is  to  remap  the  pressure  response  in  the 
wavenumber  domain  at  each  range  onto  a  new  coordinate  axis  representing  the 
moving  receiver.  This  mapping  can  be  implemented  by  a  bulk  shift  and  fine 
interpolation  of  the  complex  wavenumber  spectra  in  the  positive  frequency 
direction  for  a  closing  receiver,  or  in  the  negative  frequency  direction  for  an 
opening  receiver.  A  separate  interpolation  across  frequency  must  be  performed 
at  each  wavenumber,  since  the  amount  of  shifting  will  change  in  a  manner  similar 
to  the  plots  presented  in  Fig.  2a. 

To  describe  the  interpolation  procedure,  let  us  first  define  a  new  frequency 
axis  that  accounts  for  both  source  and  receiver  motion 

/,"=/r-^V"+(<'-l)A/"  ,  i  =  ,  (48) 

where 

A/'"=(/"n.„-/"™.,)/A'/  •  (49) 

Here,  the  maximum  and  minimum  frequencies  have  been  adjusted  to  account  for 
the  max  and  min  Doppler  shift  associated  with  receiver  motion.  The  Doppler 
shifts  for  each  wavenumber,  and  their  extrema,  and  ,  are  defined 
as 
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(50) 


AF„  =/^^cos(^„ 

Co 

^mm=min(AFj  ,  (51) 

n 

^max  =  max(AF„)  \fn  =  \...N^  ,  (52) 

n 


where  Vr  is  the  speed  of  the  receiver,  and  is  its  direction  relative  to  grazing. 

The  maximum  or  minimum  frequency  shifts  are  used  to  modify  the  edge 
values  of  the  new  frequency  axis  depending  on  the  direction  of  receiver  motion 
(i.e.,  positive  away  from  the  source  or  negative  towards  the  source),  as  indicated 
by 


j/max+AiVm  ,  <  0 

1  /„ax  ’  Otherwise 


(53) 


|/min+^max  ,  ' 

1  /min  .  Otherwise 


(54) 


where  /,„i„and  /^ax^ere  defined  previously  in  Eqs.  (43)  and  (44)  for  the  zero 
receiver  speed  case,  and  and  defined  by  Eqs.  (51)  and  (52). 


Given  the  new  frequency  axis,  a  new  PE  field  in  the  wavenumber  domain 
can  be  generated  from  the  following  interpolation  across  frequency. 


M 

w"  (i'Ar,  KJ,")  =  ^¥'  (yAr,  k„ ,  ^  j  •  /?„ ,  (m) 

m=l 


(55) 


where  Bn  is  an  integer  shift  in  frequency  for  each  wavenumber  kn,  and  h„j(m)  are 
the  coefficients  of  a  filter  whose  group  delay  is  equivalent  to  the  remainder 
between  the  actual  required  frequency  shift  and  the  bulk  shift  B„.  The  B„  values 
are  defined  as  B^  =  FLOOR[^F^\,  using  the  MATLAB  command  as  an  operator. 

The  coefficients  of  hn,i(m)  can  be  designed  using  procedures  similar  to  those 
applied  to  time  domain  interpolative  beamforming.  For  the  case  of  a  two 
coefficient  linear  interpolator  (used  in  our  implementation),  the  values  of  hn,i(m) 
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can  be  deterministically  defined  as 


(56a) 

(56b) 


It  is  again  noted  that  the  interpolation  is  performed  separately  across  the 
frequency  axis  for  each  discrete  wavenumber  coordinate  at  each  range  cell.  It  is 
also  noted  that  the  values  y/"{jl^r,k^,f")  represent  the  final  output  field,  but  the 

un-interpolated  field  y/'{jAr,k^,f.')  is  propagated  to  the  next  range  within  the 
marching  loop. 
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III.  M3PE  IMPLEMENTATION  DETAILS 


To  complete  our  description  of  the  propagation  model,  we  also  need  to 
discuss  approaches  to  handling  a  variety  of  implementation  details.  These 
include  a)  taper  functions  to  handle  finite  computational  depths,  b)  transition 
functions  for  sound  speed  and  density  discontinuities  at  the  interface  between 
water  and  the  bottom,  c)  incorporation  of  attenuation,  d)  options  for  bottom 
roughness,  e)  rules  of  thumb  for  range  and  depth  mesh  size  resolutions,  and  f) 
known  limitations  of  the  current  M3PE  implementation.  Many  of  the  limitations 
stated  here  have  already  been  implemented  in  Smith’s  MMPE  model,  so 
incorporating  them  into  the  MATLAB  version  is  fertile  territory  for  future 
upgrades. 

A.  TAPER  FUNCTIONS 

In  the  M3PE  implementation,  the  computational  depth,  Zmax,  is  nominally 
twice  the  maximum  water  depth,  Zw.  To  minimize  the  effects  of  energy 
propagating  at  very  deep  sub-bottom  depths,  or  from  very  steep  vertical  angles, 
a  filter  is  applied  to  both  the  wavenumber  (KE)  and  depth  (PE)  propagators.  A 
depth  filter  is  applied  to  ensure  that  the  field  amplitude  approaches  zero  at  the 
maximum  computational  depth  (below  the  bottom  of  the  water  column), 
consistent  with  the  far-field  radiation  boundary  condition,  p{z^co)^0.  Smith^° 
points  out  that  the  filter  must  be  a  smooth  function  to  avoid  generating  higher 
order  spectral  terms  during  Fourier  transformation  from  the  abruptness  of  the 
filter  itself.  According  to  Smith,^°  it  is  desired  that  the  filter  begin  its  taper 
approximately  1/3  of  the  distance  from  the  Zmax-  A  function  that  creates  the 
desired  affect  is  defined  as 

G(z)=i(l-cos(n,))+i  ,  (57) 

with 
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Q,  =  min 


(58) 


i  =  • 


2  ’  2 


As  defined  previously,  is  the  size  of  the  Discrete  Fourier  Transform  (DFT).  In 
the  model  implementation,  the  depth  propagator  U{z)  (see  Eq.  (17))  is  multiplied 
by  G(z)  at  a  one-time  initialization  step.  The  filter  is  then  repeatedly  applied  at 
each  range  step,  resulting  in  severe  attenuation  of  the  deepest  depths  after 
several  updates  of  the  marching  loop.  A  plot  of  the  depth  taper  function  is  shown 
in  Fig.  3a. 

To  restrict  propagation  at  very  steep  angles  (where  the  PE  approximation 
is  not  valid)  a  separate  function  is  applied  to  the  wavenumber  propagator  7op(A:J 
(see  Eq.  (20))  of  the  form 


^(«)  =  ^(l  +  cos(0j) 

where 


0„  =7rx< 


0 


<  sin  ( >  sin  ( ) 

Kq  Kq 

—  >  sin  (0  ) 

Kq 

Kq 


(59) 


(60) 


with  the  control  parameters  nominally  defined  as  =  80“  and  0^,  =  90“ . 

These  parameters  result  in  a  very  sharp  transition  function  as  illustrated  in 
Fig.  3b.  This,  according  to  Smith, provided  the  best  results  during  benchmark 
testing.  A  smaller  value  for  0^^  would  result  in  a  smoother  response  but  would 
not  excite  the  higher  angles  of  propagation  as  well. 
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Wavenumber  taper  Auction 


ampitiKle 


(a)  (b) 

Figure  3.  Examples  of  a  depth  taper  (a)  and  wavenumber  taper  (b)  for  the  case 

f=100  Hz,  Co=1500  m/s,  and  D=200  m 


B.  BOTTOM  TRANSITION  FUNCTIONS 

When  using  a  Discrete  Fourier  Transform  to  perform  a  spectral 
decomposition,  it  is  standard  practice  to  smooth  any  sharp  amplitude  transitions 
in  the  data  to  avoid  creating  harmonic  components  representing  the  true  Fourier 
Series  of  any  finite  width  envelope.  In  the  case  of  the  split-step  marching 
algorithm,  these  “harmonics”  would  be  observed  as  sidelobes  along  the  depth 
axis.  In  the  parametric  data  representing  an  acoustic  environment,  it  is  likely  that 
a  sharp  transition  will  occur  in  the  sound  speed  and  density  values  at  the  water- 
bottom  interface.  In  Smith’s  MMPE  and  our  M3PE  implementation,  a  separate 
transition  function  is  applied  to  the  sound  speed  and  the  density.  For  the  speed 
of  sound,  the  transition  function  is  defined  as 


hM) 


1  +  e  ' 

V  y 


where 


and 


L„  =  max 


— ,Az 

vlO 


(61) 


(62) 


(63) 
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Note  that  the  maximum  function  in  Eq.  (63)  enforces  the  constraint  that  4>Az. 
Hc(CJ  has  the  following  properties: 


hAC)  = 


<  0.5 
^1 


z<z, 
z  =  z, 
z  >  z, 


(64) 


The  transition  function  in  Eq.  (61)  is  applied  to  the  sound  speed  at  each  range 
about  the  bottom  depth  according  to 


c{z,r)  =  c{z,r)  +  {c^-c„)-HXz-Zb)  ■ 


(65) 


A  different  function,  having  similar  properties  to  He,  is  applied  to  the 
density  data.  This  transition  function  is  a  cubic  spline  defined  as 


1 

/  V 
1  +  ^ 

V 


2  3 


kAj 


1-2 

3 


i-i 

V  i./ 


-L  L 

2  2 


(66) 


2 


where 


Lp  =  max(2/l,2Az) 


(67) 


H^Q  is  applied  to  the  density  parameters  p(z,r)  using  Eq.  (65)  with  the 
density  and  density  transition  function  substituted  for  the  appropriate  sound 
speed  related  parameters.  According  to  Smith, Hp  was  chosen  to  be  a 
continuously  defined  function  to  guarantee  the  existence  of  its  second  derivative. 
Recall  that  the  propagator  U2  in  Eq.  (29)  is  proportional  to  i-e.  the 

second  derivative  of  the  transition  function.  Smith  defined  this  function  as 
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»:(()= 


L- 


1  +  ^ 

V  ^pj 
f 


-4 


\^p  J 

x-i- 

V  J 


-L  L 

_^<r<^ 

2  2 


-L 

— ^<C<L 
1  P 


(68) 


Examples  of  Hp  and  //"  are  plotted  in  Fig.  4  for  a  120  m  waveguide  at  a 
frequency  of  100  Hz. 


Sound  speed  transition  function  Density  transition  function  2nd  Derivative  Dens  Transition  fen 


Hc(z-z,) 


H,"(z-z>,) 


X  10 


*3 


(a)  (b)  (c) 

Figure  4.  Example  water-bottom  transition  functions  for  (a)  sound  speed  and  (b) 
density  when  =  15  m,  and  (c)  2'^^  derivative  of  density  function. 


C.  ATTENUATION 

Attenuation  is  implemented  by  adding  an  exponential  damping  term  to  the 
depth  propagator  function  used  in  the  PE  marching  algorithm,  specifically 


-ik,,—Unp(f+Ar,z)-—a(r+Ar,z) 

i//{r  +  Ar,z)  =  e  ^  ^  xIDFTle 


ikaArtop(k^] 


xDFT 


-ik„^Uop(r,z)~air,z) 

e  ^  ^  xy/{r,z) 


,  (69) 
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where  a(r,z)  =  a(r,z)/8.686,  in  units  of  nepers,  and  a  is  the  attenuation  in  units  of 
dB/m. 


The  attenuation  is  defined  separately  for  the  water  column  versus  the 
bottom  layer.  In  the  water,  our  implementation  employs  the  Fisher-Simmons 
model  as  defined  by  Kinsler,  et  al.,^^  given  by 


a(r,z) 


A 


B 


■  +  c 


f 


1000 


in  dB/m 


(70) 


where 


^  =  0.083 


V35y 


c  = 


(4.9x10 


T  T 

/i  =  780^29  ,  =  42000e^  ,  and 


(71) 

(72) 

(73) 

(74) 

(75) 


Although  subscripts  were  omitted  for  convenience,  the  environmental 
parameters  T  (temperature  in  °C),  S  (salinity  in  ppm) ,  and  pH  (acidity  coefficient, 
dimensionless)  are  allowed  to  be  range-dependent.  However,  in  the  current 
implementation  of  M3PE,  only  a  single  value  of  T,  S,  and  pH  is  used  to  compute 
attenuation  in  the  entire  water  column  at  each  range  cell.  In  the  bottom  layer,  the 
attenuation  is  simply  defined  as  a  range  dependent  constant,  a\r,z),  in  terms  of 
dB/A..  This  is  converted  to  dB/m  according  to 

a(r,z)=  a' (r,z)—^  .  (76) 

c{r,z) 
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D.  BOTTOM  ROUGHNESS 

The  M3PE  implementation  allows  the  mean  ocean  depth  versus  range  to 
be  modified  by  a  range-dependent  roughness  function.  Several  options 
implemented  in  M3PE  include  a  zero  mean  Gaussian  and  uniform  random 
perturbation  with  a  user  defined  RMS  height.  The  better  option  is  based  upon  the 
original  work  of  Fox  and  Hayes^^  where  the  seafloor  roughness  is  modeled  with 
an  isotropic,  zero  mean,  magnitude  spectrum  that  follows  a  power  law.  Our 
implementation  is  based  upon  the  description  given  by  Smith^"^  as  outlined  below. 

First  define  the  envelope  of  a  roughness  spectrum  having  the  form 


fV(n)  = 


r?  —  1 ,2,3,  . . . ,  A/f 


(77) 


where  Lcor  is  a  correlation  length,  /c„  is  the  horizontal  wavenumber.  Nr  is  the 
number  of  range  cells,  and  p  is  a  power  constant  that  is  usually  set  to  the  value 
of  3.5.  The  spectral  envelope  is  converted  into  a  bottom  depth  perturbation  by 
using  Eq.  (77)  to  factor  a  complex  exponential  with  a  Gaussian  amplitude  and 
uniform  phase.  To  illustrate,  let 

5(n)  =  1,2,3,  ...,/V,  ,  (78) 


where  A„  is  a  standardized  Gaussian  random  variable  with  zero  mean  and  unit 
standard  deviation  (i.e.,  A/(0,1)),  6’„  is  a  uniform  random  variable  between  the 
interval  [0,1],  and  Ak  is  the  wavenumber  cell  resolution  equal  to  2nlRmax.  Rmax  is 
the  length  of  the  range  interval.  Then  the  spectrum  function  is  transformed  into  a 
range  dependent  vector  by 


^\j)  =  ^S{n)e 


,j=  1,2,3,  ...,  Nr 


n=l 


(79) 


Next,  Eq.  (79)  is  normalized  and  scaled  to  result  in  the  final  perturbation 

vector. 
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^(7)  = 


1 

2 


,j  =  1,2,3,  ...,  Nr 


(80) 


i'  N  \ 

f^W^n) 


V«=1 


where  a  is  a  user  defined  scale  factor. 

The  perturbation  vector  specified  by  Eq.  80  is  the  same  length  as  the 
range  grid.  This  vector  is  added  to  the  mean  bottom  depth  at  each  discrete  range 
to  result  in  a  “roughened”  bottom.  Thus  the  range  resolution  imparts  a  lower  limit 
on  the  granularity  of  roughness  features.  An  example  of  a  roughness  vector 
computed  for  a  point  set  of  parameters  is  shown  in  Fig.  5. 


Example  bottom  roughness 


Figure  5.  Example  Bottom  Roughness 


E.  RANGE  AND  DEPTH  MESH  SIZES 

In  reference  [20],  Smith  evaluated  the  convergence  properties  of  various 
depth  and  range  mesh  sizes  for  various  benchmark  problems.  It  was  found  that 
optimal  performance  was  obtained  with  the  approximate  scales 


Az«— ,  (81) 

10 

.  (82) 


where  X  is  the  signal  wavelength. 
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The  requirement  for  depth  resolution  is  dominated  by  the  sound  speed 
mixing  length  definition,  which  is  max(Az,A,/10).  Any  depth-mesh  size  greater 
than  yiO,  will  further  increase  the  transition  distance  of  the  sound  speed  at  the 
bottom,  thereby  causing  inaccuracies  in  the  model  solution.  The  range  step  size 
appears  to  be  determined  by  the  “optimal”  amount  of  phase  information  that  is 
propagated  at  each  step  in  the  marching  algorithm.  Too  large  a  step  size  causes 
errors  due  to  stationarity  limits,  while  too  small  a  step  size  introduces  errors 
associated  with  numerical  accuracy^°.  Chapter  IV  will  show  examples  of 
performance  on  benchmark  environments  for  various  range  step  sizes. 


F.  LIMITATIONS  OF  M3PE  IMPLEMENTATION 

It  is  important  to  note  that  the  M3PE  code  does  not  implement  all  the 
features  in  Smith’s  original  MMPE  code,  as  well  as  a  variety  of  other  useful 
capabilities  that  could  be  included  in  future  upgrades.  This  section  contains  a 
short  description  of  these  limitations. 

Currently,  the  M3PE  code  only  supports  the  modeling  of  a  sinusoidal 
source  signal.  Signals  with  other  spectral  shapes  could  easily  be  incorporated  by 
upgrading  the  source  definition  function  to  create  the  broadband  signal  and 
appropriately  adjusting  each  frequency  component  to  reflect  the  required  Doppler 
at  each  propagating  wavenumber.  The  current  MMPE  model  supports  a 
broadband  spectrum  with  the  shape  of  a  Hanning  envelope. 

Next,  the  bottom  fluid  model  in  M3PE  currently  does  not  represent  the 
effects  of  shear.  As  noted  by  Smith^°,  this  can  be  accomplished  by  creating  an 
effective  density  and  attenuation  that  includes  the  shear  sound  speed,  Cs,  and 
shear  attenuation,  as,  parameters.  Another  feature  not  supported  by  M3PE  is  any 
treatment  of  radial  coupling.  Currently,  only  a  single  radial  (range  versus  depth) 
response  is  generated. 

Also,  no  accommodation  is  provided  for  more  than  one  bottom  layer. 
Sometimes,  it  is  desirable  to  represent  a  “slow”  sediment  layer  above  a  hard  (i.e.. 
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fast)  sub-bottom  layer.  This  can  be  implemented  by  applying  the  discontinuity 
functions  for  sound  speed  and  density  to  the  transition  regions  at  the 
water/sediment  depths  as  well  as  at  the  sediment/sub-bottom  interface. 


Finally,  no  allowance  has  been  made  for  surface  roughness.  One 
technique,  first  introduced  by  Dozier^^,  is  based  upon  a  conformal  mapping 
algorithm  that  transforms  the  depth  variant  variables  (e.g.,  c{r,z),  p{r,z))  normally 
affected  by  the  variable  surface  height,  into  a  pseudo-space  that  has  a  flat 
surface.  In  the  SSF  method,  the  parameters  in  the  pseudo-space  are  used  to 
propagate  the  solution  for  one  range  step,  and  then  inverted  back  to  physical 
space  for  the  next  iteration.  The  result  is  a  propagated  field  resulting  from  a 
single  realization  of  an  irregular  surface.  Another  technique  presented  by  Tappert 
and  Nghiem-Phu^^  uses  the  method  of  images  to  solve  the  PE  with  an  rough 
surface.  Either  of  these  methods  would  be  a  prime  candidate  for  future  upgrades. 
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IV.  M3PE  APPLIED  TO  BENCHMARK  CASES 


Prior  to  using  the  M3PE  code  to  investigate  the  frequency  spread  of 
sinusoidal  signals  in  the  presence  of  motion,  the  performance  of  the  M3PE  model 
was  validated  against  a  set  of  known  test  cases.  A  well  documented  and 
sufficiently  challenging  set  of  benchmark  environments  is  conveniently  available 
from  the  website  associated  with  the  Shallow  Water  Acoustic  Modeling  workshop 
held  in  1999  (a.k.a.  SWAM’99)^^.  Several  of  the  SWAM’99  environments  were 
processed  and  compared  to  results  generated  by  Smith  using  the  original  MMPE 
code,  and  by  Mikhin  who  used  an  energy  conserving  implicit  finite  difference 
(IFD)  PE  model^^.  The  IFD  models  are  generally  considered  to  be  more  accurate 
than  SSF  algorithms  because  of  their  higher  order  treatment  of  the  bottom 
interface  boundary  condition. Our  intention  is  to  demonstrate  close  agreement 
between  the  M3PE  and  MMPE  implementations,  as  well  as  to  make  general 
comments  regarding  how  well  either  SSF  implementation  approaches  the 
performance  of  an  IFD  method.  Since  our  intention  here  is  to  validate  the 
execution  of  the  core  SSF  PE  implementation,  both  the  source  and  the  receiver 
have  zero  speed  in  all  the  examples  in  this  chapter. 

A.  TRANSMISSION  LOSS  FOR  FLATa 

The  first  environment,  called  FLATa,  is  characterized  by  a  flat  bottom  and 
a  constant  sound  speed  (Cw)  and  density  (pw)  in  the  water  column.  A  ray  trace 
illustrating  iso-speed  propagation  in  the  water  column  for  the  FLATa  environment 
is  shown  in  Fig.  6,  with  a  source  depth  of  30  m.  The  bottom  parameters  for  this 
case  are  quite  range-dependent  along  a  20  km  radial  distance.  Illustrations  of  the 
range-dependent  bottom  parameters  including  bottom  sound  speed  (Cb),  gradient 

of  bottom  sound  speed  bottom  density  (pb)  are  plotted  in  Fig.  7.  The 

attenuation  is  zero  in  the  water  column  (aw),  and  is  a  constant  0.1  dB/A,  in  the 
bottom  (ab).  The  shear  sound  speed  (Cs)  and  attenuation  (as)  were  assumed  to 
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be  zero.  A  color  coded  representation  of  sound  speed,  density  and  attenuation 
versus  range  and  depth  is  provided  in  Fig.  8.  It  can  be  observed  that  the 
attenuation  values  follow  the  shape  of  the  bottom  sound  speed  to  maintain  a 
constant  loss  per  wavelength.  A  summary  of  all  parameters  is  given  in  Table  1 . 


Figure  6.  Ray  trace  in  water  column  for  the  FLATa  environment,  Zs  =  30  m 


Figure  7.  Bottom  parameters  vs  range  for  the  FLATa  environment 
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Figure  8.  (a)  Sound  speed,  (b)  density,  and  (c)  attenuation  vs  range  and  depth 

for  the  FLATa  environment 
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An  example  transmission  loss  (TL)  from  the  M3PE  model  computed  for 
the  FLATa  environment  at  250  FIz  is  shown  in  Fig.  9.  The  operating  parameters 
such  as  DFT  size  (A/z),  number  of  range  cells  (Nr),  computational  depth  (Zmax), 
and  source  depth  (Zs)  are  listed  beside  the  range  vs  depth  plot.  At  250  Hz,  the 
wavelength  {X)  is  about  6  m.  Thus  the  range  resolution  (5.0  m)  and  depth 
resolution  (0.195  m)  are  close  to  the  settings  recommended  by  Smith^°  as 
defined  in  Eqs.  (81)  and  (82). 
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Transmission  loss  (dB),  f  =  250  Hz.  z  =  30  m,  v  =0  Kts.  v  =  0  Kts 
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Nz  =  2048 
Nr  =  4000 
Az  =  0.1953  m 
Ar  =  5.0  m 
Zb  =  100  m 
Zmax  =  200  m 
f=  250  Hz 
Zs  =  30  m 


Figure  9.  Transmission  loss  vs  range  and  depth  at  250  Hz  for  the  FLATa 

environment 


A  comparison  of  the  TL  generated  by  M3PE,  MMPE,  and  the  IFD  models 
is  plotted  in  Fig.  10  for  a  250  Hz  source  frequency,  a  source  depth  of  30  m,  and  a 
receiver  depth  of  35  m.  The  top  plot  in  Fig.  10  is  TL  for  the  entire  20  km  range. 
The  second  plot  is  an  expansion  of  the  first  5  km  while  the  third  plot  is  an 
expansion  of  the  last  5  km.  It  can  be  observed  in  Fig.  10  that  the  M3PE  (blue) 
and  MMPE  (red)  results  are  very  close  in  the  first  5  km  and  reasonably  close  at 
the  latter  ranges.  Also,  all  three  models  follow  quite  closely  out  to  about  2700  m, 
but  somewhat  diverge  towards  the  latter  ranges.  Still  the  top  plot  in  Fig.  10 
illustrates  that  the  general  trend  of  the  TL  response  is  very  similar  even  at  the 
longer  ranges. 
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Zs  =  30  m,  Zr  =  35  m,  M3PE  ^  Az  =  0.1953,  Ar  =  5.0  m,  MMPE  Az  =  0.1953,  Ar  =  5.0  m 


Figure  10.  Comparison  of  transmission  loss  at  250  Hz  between  MMPE,  Mikhin’s 

IFD,  and  M3PE  PE  models  for  FLATa  environment 


Figure  1 1  presents  a  TL  prediction  from  the  M3PE  model  at  25  Hz.  Here 
Nz  =  256  and  Nr  =  400,  consistent  with  Smith’s  guidance  regarding  grid  sizes.  It 
can  be  observed  that  significantly  more  bottom  penetration  is  predicted  at  this 
frequency  than  observed  at  250  Hz,  as  expected.  A  comparison  between  TL 
results  generated  by  the  three  different  models  at  25  Hz  is  plotted  in  Fig.  12.  In 
this  case,  the  M3PE  and  MMPE  results  are  quite  similar  all  the  way  to  the 
maximum  range.  However,  both  SSF  models  differ  significantly  from  the  IFD 
approach  after  the  first  2500  m.  Smith^^  explains  that  this  dissimilarity  is  related 
to  the  bottom  mixing  functions  applied  to  the  sound  speed  and  density,  which  are 
proportional  to  either  the  wavelength  or  the  depth  mesh  size  (see  Eqs.  (63)  and 
(67)).  At  low  frequencies  the  mixing  functions  can  be  quite  wide  (in  depth)  and 
cause  fairly  large  changes  in  the  actual  cjct,  and  pw/pb  parameters  near  the 
bottom.  For  shallow  waveguides  this  can  represent  a  significant  deviation.  This  is 
a  known  disadvantage  of  SSF  approaches,  and  thus  limits  their  effectiveness 
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when  the  ocean  depth  is  a  very  small  number  of  wavelengths. 


Nz  =  256 
Nr  =  400 

Az  =  1 .56  m 
Ar  =  50.0  m 
Zb  = 100  m 
Zmax  =  200  m 
f  =  25  Hz 
Zs  =  30  m 
Zr  =  35  m 


Figure  1 1 .  Transmission  loss  vs  range  and  depth  at  25  Hz  for  the  FLATa 

environment 


Figure  12.  Comparison  of  transmission  loss  at  25  Hz  between  MMPE,  Mikhin’s 
IFD,  and  M3PE  PE  models  for  the  FLATa  environment 
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Additionally,  comparisons  were  made  with  the  FLATa  environment  at  1000 
Hz.  Figure  13  illustrates  the  2-D  transmission  loss  and  Fig.  14  plots  the  results  at 
Zs  =  30  m  and  Zr  =  35  m  produced  by  the  three  different  models.  Consistent  with 
the  250  Hz  case,  the  TL  produced  by  M3PE  and  MMPE  at  1000  Hz  are  quite 
similar,  and  both  follow  the  general  trends  generated  by  the  IFD  algorithm. 

Transmission  loss  (dB),  f  =  1000  Hz,  =  30  m,  =  0  Kts,  =  0  Kts 
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„  Nz  =  2048 
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_eo  Zs  =  30  m 
Zr  =  35  m 
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Figure  13.  Transmission  loss  vs  range  and  depth  at  1000  Hz  for  the  FLATa 

environment 
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Zs  =  30  m,  Zr  =  35  m,  M3PE  Az  =  0.195,  Ar  =  2.5  m,  MMPE  Az  =  0.193,  Ar  =  10.0  m 

Figure  14.  Comparison  of  transmission  loss  at  1000  Hz  between  MMPE,  Mikhin’s 

IFD,  and  M3PE  PE  models  for  the  FLATa  environment 

B.  CONVERGENCE  DEMONSTRATED  ON  FLATa 

The  FLATa  environment  was  also  used  to  evaluate  the  convergence 
properties  of  the  M3PE  code.  Here  we  are  specifically  interested  in  the  relative 
performance  as  a  function  of  the  number  of  range  cells  (Nr)  which  determines  the 

range  resolution  Ar  =  ^^.  Smith  showed  that  an  SSF  PE  model  can  produce 

significantly  different  results  for  various  resolutions.  This  effect  is  of  interest  to 
this  study  because  processor  and  memory  resources  can  be  highly  stretched 
when  the  number  of  range  cells  is  large.  This  is  especially  important  when  the 
model  is  evaluated  across  many  frequencies,  as  will  be  required  when  we  model 
the  effects  of  source  and  receiver  motion.  The  following  examples  demonstrate 
that  the  convergence  properties  of  M3PE  are  similar  to  the  results  obtained  with 
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the  MMPE  implementation.  We  also  show  that  the  guidance  provided  by  Smith^°, 
and  restated  in  Eqs.  (81)  and  (82),  is  very  consistent  with  the  results  presented 
here. 

A  comparison  of  TL  vs  range,  in  the  FLATa  environment,  using  various 
range  resolutions  (Ar  =  2.5  m,  Ar  =  5.0  m,  Ar  =  10  m,  Ar  =  20  m)  is  provided  in 
Fig.  15.  Again,  the  source  depth  is  30  m  and  the  receiver  depth  is  35  m.  The 
frequency  is  250  Hz,  and  consequently  the  wavelength  is  6  m.  Here  we  see  that 
a  resolution  of  at  least  5  m  (A^,.  =  2000)  is  required  for  convergence  at  the  longer 
ranges.  This  is  quite  apparent  where  the  TL  with  the  lowest  resolution  can  be 
several  dB  different  than  the  TL  with  the  “optimal”  resolution.  It  is  interesting  that 
there  is  not  much  degradation  with  lower  resolution  inside  5  km. 

Figure  16  contains  a  second  test  case  with  the  FLATa  environment  at 
1000  Hz.  The  transmission  loss  is  again  plotted  for  the  same  four  resolutions. 
Here  we  see  that  at  this  higher  frequency,  with  X  =  1.5  m,  there  is  quite  a  bit  of 
difference  between  the  TL  response  using  the  four  different  range  resolutions. 
Presumably  the  more  accurate  result  is  when  the  range  resolution  closely 
matches  the  wavelength  (i.e.  Nr  =  8000).  Similar  to  the  250  Hz  case,  it  appears 
there  is  a  loss  of  energy  when  the  resolution  is  less  than  required.  However, 
there  isn’t  a  dramatic  difference  until  resolution  is  more  than  2-3  times  the 
wavelength. 
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Figure  15.  Comparison  of  M3PE  results  at  250  Hz  for  the  FLATa  environment,  Ar 
=  2.50  m  (magenta),  Ar  =  5.0  m  (blue),  Ar  =  10.0  m  (red),  Ar  =  20.0  m 

(green) 


Figure  16.  Comparison  of  M3PE  results  at  1000  Hz  for  the  FLATa  environment, 
Ar  =  2.50  m  (magenta),  Ar  =  5.0  m  (blue),  Ar  =  10.0  m  (red),  Ar  =  20.0 

m  (green) 
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C.  RECIPROCITY  DEMONSTRATED  ON  FLATa 

The  FLATa  environment  was  also  used  to  evaluate  the  reciprocity 
characteristics  of  the  M3PE  model.  Here  we  evaluated  the  TL  beginning  with  a 
source  at  30  m  and  a  receiver  20  km  away  at  a  depth  of  35  m.  The  frequency  is 
250  Hz,  the  depth  resolution  is  0.1953  m,  and  the  range  resolution  is  5  m.  Here, 
the  TL  is  propagated  from  “left  to  right”  across  the  environment  illustrated  in  Fig. 
8.  Conversely,  another  TL  is  computed  starting  at  a  35  m  depth  and  propagating 
from  right  to  left  across  the  environment  to  be  received  at  a  depth  of  30  m.  The 
TL  (level  and  phase)  associated  with  the  last  km  of  distance  is  plotted  on  the 
same  axis  in  Fig.  17.  The  level  and  phase  at  20  km  resulting  from  propagation 
across  reciprocal  paths  is  practically  identical.  The  performance  demonstrated 
here  is  fairly  similar  to  the  results  generated  by  Smith^°. 


Forward  ^  r  =  0  m,  Zs  =  30  m,  z,  =  35  m,  Reciprocal  ^  r  =  20  km,  Zs  =  35  m,  z,  =  30  m 


Figure  17.  Results  of  a  reciprocity  test  for  the  FLATa  environment,  f  =  250  Hz. 

Magnitude  of  TL  (upper)  and  phase  (lower) 

D.  TRANSMISSION  LOSS  FOR  DOWNa 

The  second  benchmark  is  called  DOWNa.  This  environment  is 


39 


characterized  by  an  increasing  ocean  depth  versus  range,  and  a  constant  sound 
speed  (Cw)  and  density  (pw)  in  the  water  column.  A  ray  trace  illustrating  iso-speed 
propagation  in  the  water  column  in  a  range-dependent  waveguide  is  shown  in 
Fig.  18  with  a  source  depth  equal  to  30  m.  Except  for  ocean  depth,  the  other 
bottom  parameters  such  as  bottom  sound  speed  (Cb),  gradient  of  bottom  sound 
/ 

speed  (  ),  and  bottom  density  (pb)  are  constant  versus  range.  Also,  the 

attenuation  is  zero  in  the  water  column  (aw),  and  is  a  constant  0.1  dB/?^  in  the 
bottom  (ab).  Again  the  shear  sound  speed  (Cs)  and  attenuation  (as)  are  zero.  A 
color  coded  representation  of  sound  speed,  density  and  attenuation  versus  range 
and  depth  for  this  environment  is  provided  in  Fig.  19.  As  before,  the  attenuation 
values  follow  the  shape  of  the  bottom  sound  speed  to  maintain  a  constant  loss 
per  wavelength.  A  summary  of  all  parameters  for  the  DOWNa  environment  is 
given  in  Table  2. 


Figure  18.  Ray  trace  in  water  column  for  the  DOWNa  environment,  Zg  =  30  m 
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Figure  19.  Sound  speed  (a),  density  (b),  and  attenuation  (c)  vs  range  and  depth 

for  the  DOWNa  environment 


Table  2.  Environmental  parameters  for  the  DOWNa  environment 
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An  example  TL  from  the  M3PE  model  against  the  DOWNa  environment  at 
250  Hz  is  shown  in  Fig.  20.  The  operating  parameters  Nz  and  Nr  are  listed  beside 
the  TL  plot.  In  this  case,  the  range  resolution  meets  the  criterion  of  Eq.  (82),  but 
the  depth  resolution  is  approximately  double  the  recommended  guidance  from 
Eq.  (81).  A  comparison  of  the  TL  at  250  Hz  generated  by  M3PE,  MMPE,  and  the 
IFD  models  is  plotted  in  Fig.  21.  It  can  again  be  observed  that  the  M3PE  (blue) 
and  MMPE  (red)  results  are  very  close  in  the  first  5  km  and  reasonably  close  at 
the  latter  ranges.  Also,  the  general  trend  of  the  TL  response  from  the  two  SSF 
methods  follows  the  trend  of  the  IFD  approach  out  to  12  km,  except  for  a  few 
anomalous  ranges. 
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f  =  250  Hz 
Zs  =  30  m 
Zr  =  35  m 


Figure  20.  Transmission  loss  vs  range  ancJ  cJepth  at  250  Hz  for  the  DOWNa 

environment 


Figure  21 .  Comparison  of  transmission  loss  at  250  Hz  between  MMPE,  Mikhin’s 
IFD,  and  M3PE  PE  models  for  the  DOWNa  environment 
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E.  TRANSMISSION  LOSS  FOR  IWa 

The  final  benchmark  is  characterized  by  a  constant  ocean  depth  and 
density  (pw)  in  the  water  column,  but  a  highly  variable  water  sound  speed  {c„). 
Since  the  range-  and  depth-dependent  sound  speed  profile  is  intended  to 
represent  the  possible  effects  from  internal  waves,  the  environment  is  labeled 
IWa.  A  ray  trace  illustrating  propagation  in  the  water  column  is  shown  in  Fig.  22 
(again  =  30  m).  All  the  other  environmental  parameters,  such  as  bottom  sound 

dc  / 

speed  {cb),  gradient  of  bottom  sound  speed  (  ),  and  bottom  density  (pb)  are 

constant  versus  range.  Again,  the  attenuation  is  zero  in  the  water  column  (aw) 
and  0.1  dB/?^  in  the  bottom  (ab).  The  shear  sound  speed  (c^)  and  attenuation  (as) 
are  also  zero.  A  color  coded  representation  of  sound  speed,  density  and 
attenuation  versus  range  and  depth  for  this  environment  is  illustrated  in  Fig.  23.  A 
summary  of  all  parameters  for  the  IWa  environment  is  given  in  Table  2. 


Ray  Trace  Plot,  6  rays,  start  angles  [-5.0,5  0]  degs.  zs  =  30  m 


Figure  22.  Ray  trace  in  water  column  for  the  IWa  environment,  Zs  =  30  m 
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sound  speed  (m/s)  in  water 


Figure  23.  SouncJ  speecJ  vs  range  ancJ  cJepth  for  the  IWa  environment 
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The  ray  trace  in  Fig.  22  shows  a  very  complicated  propagation 
environment  with  multiple  bottom  bounce  propagation  without  many  surface 
interactions  due  to  the  high  sound  speeds  near  the  surface.  Flowever,  a  weak 
surface  duct  is  also  supported  for  small  range  sections.  An  example  TL  versus 
range  and  depth  plot  is  provided  in  Fig.  24. 
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Figure  24.  Transmission  loss  vs  range  and  depth  at  250  Hz  for  the  IWa 

environment 


A  comparison  between  the  TL  generated  by  the  M3PE,  MMPE  and  the 
IFD  model  for  the  IWa  environment  is  shown  in  Fig.  25.  Again  Zs  is  30  m  and  Zr  is 
35  m.  It  can  be  observed  that,  in  this  case,  all  three  models  produce  similar 
results  out  to  about  4000  m.  Also,  the  dissimilarities  at  the  longer  ranges  are 
fairly  modest.  However,  the  complexity  of  the  sound  speed  makes  it  difficult  to 
draw  any  conclusions  regarding  the  closeness  of  detailed  fluctuations  at  the 
longer  ranges.  It  suffices  to  say  that  the  general  trend  of  all  three  models  across 
the  20  km  propagation  distance  is  fairly  similar. 
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-TL  (dB)  -TL  (dB)  -TL  (dB) 


M3PE  ^  Az  =  0.391 ,  Ar  =  5.0  m,  MMPE  ^  Az  =  0.391 ,  Ar  =  10.0  m 


Figure  25.  Comparison  of  transmission  loss  at  250  Hz  between  MMPE,  Mikhin’s 

IFD,  and  M3PE  PE  models  for  the  IWa  environment 


46 


V.  BANDWIDTH  ESTIMATION  AND  STATISTICAL 
CHARACTERIZATION 


The  linear  magnitude  of  the  transmission  loss  distributed  across  range, 
depth  and  frequency  associated  with  a  moving  source  and  receiver  serves  as  the 
domain  from  which  the  received  signal  bandwidth  is  estimated.  The  bandwidth 
parameter  is  designed  to  represent  the  spreading  of  the  TL  across  frequency 
with  a  single  parameter.  There  are  several  accepted  methods  for  defining 
bandwidth.  For  example  the  effective  bandwidth  equation  used  to  characterize 
the  transfer  function  of  digital  filters  is  a  common  method.^®  The  formulation  for 
this  estimator  is  given  as 
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Here  pi  is  the  magnitude  of  the  evaluated  magnitude  response  across  (in  our 
case)  frequency  and  the  constant  Af\s  the  discrete  separation  of  each  frequency 
bin.  Other  definitions  of  the  bandwidth  parameter  are  also  used  in  the  literature, 
such  as  treating  the  frequency  dependent  pressure  magnitude  as  a  probability 
density  function  (PDF),  and  representing  the  sample  standard  deviation  as  a 
bandwidth.'^°  The  bandwidth  parameter  for  this  method  would  be 
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Still  other  methods  use  an  iterative  search  algorithm  to  find  the  frequency  limits 
that  satisfy  some  objective  function  such  as  the  rate  of  change  of  a  normalized 


where  p  refers  to  the  sample  mean,  // 
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boxcar  integration  statistic,  or  the  deflection  between  the  pressure  magnitude  at 
its  peak  and  the  lower  values  along  the  skirts. 

For  our  study,  we  choose  the  bandwidth  parameter  described  by  Eq.  (83) 
since  our  experience  has  shown  this  equation  to  give  intuitively  appealing  results 
regardless  of  the  shape  of  the  underlying  function.  For  example,  the  statistic 
equates  to  unity  when  the  spectrum  function  is  perfectly  flat.  An  example  of 
results  generated  with  this  equation  on  three  different  magnitude  responses  are 
shown  in  Fig.  26  below.  It  can  be  observed  that  the  algorithm  computes  an 
“effective”  bandwidth  representing  the  fractional  magnitude  enclosed  by  the 
magnitude  response  across  the  evaluated  area.  Thus,  the  same  bandwidth  value 
can  result  from  curves  with  remarkably  different  shapes.  However,  the  statistic 
can  provide  reasonable  parameter  estimates  even  when  the  shape  is  fairly 
complicated. 
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Figure  26.  Example  Bandwidths  Computed  for  Three  Different  Functions 

In  the  next  section,  we  compute  the  BW  parameter  from  the  frequency 
variant  transmission  loss  represented  on  a  discrete  grid  of  range  and  depth  cells 
associated  with  several  different  ocean  environments.  The  results  are  plotted  as 
an  aggregate  histogram  including  estimates  of  the  sample  mean  and  standard 
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deviation.  Additionally,  the  bandwidth  variables  are  marginalized  over  one 
independent  variable  (e.g.,  range  or  depth)  and  plotted  against  the  other.  Finally, 
scatter  plots  of  bandwidth  versus  the  normalized  linear  magnitude  of  the 
minimum  transmission  loss  across  frequency  are  generated.  The  histogram 
describes  the  overall  behavior  of  the  random  variable  across  the  entire  range- 
depth  space.  The  marginalized  plots  show  the  behavior  versus  range  or  depth, 
averaged  across  the  other  independent  variable.  Finally,  the  two-dimensional 
scatter  plots  portray  the  parameter  as  a  function  of  normalized  TL.  Any  general 
dependencies  upon  local  TL  level  are  highlighted  in  this  plot. 

The  bandwidth  statistic  is  a  measure  of  the  properties  of  the  received 
signal.  It  is  assumed  that  the  signal  being  measured  has  sufficient  signal  to  noise 
ratio  (SNR)  to  support  detection  and  accurate  parameter  estimation.  Thus  in  all 
cases  where  we  use  the  response  of  the  transmission  loss  to  estimate 
bandwidth,  we  are  assuming  that  sufficient  SNR  remains  after  the  TL  is 
accumulated  with  the  source  level  (SL),  noise  level  (NL)  and  other  components 
of  the  sonar  equation.  When  the  TL  is  so  high  to  preclude  detection  of  even  the 
loudest  known  sources  at  nominal  ranges,  the  predicted  bandwidth  parameters 
have  little  practical  value.  We  shall  see  a  case  in  the  next  section  where  this 
situation  becomes  important. 
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VI.  NUMERICAL  EXPERIMENTS  WITH  A  MOVING  SOURCE 

AND  RECEIVER 


In  this  section,  we  demonstrate  the  model  on  a  set  of  notional  ocean 
environments.  Our  intention  is  to  explore  the  way  in  which  transmission  loss 
varies  across  frequency  as  a  function  of  range  and  depth  when  both  the  source 
and  receiver  are  in  motion.  The  bandwidth  parameter  is  used  to  quantify  the 
shape  of  the  TL  vs  frequency,  including  the  dependencies  of  this  shape  on 
range,  depth,  and  TL  magnitude.  A  description  of  the  acoustic  parameters 
associated  with  the  ocean  environments  is  given  in  Section  A.  The  results 
obtained  from  the  TL  model  are  given  in  Section  B. 


A.  DESCRIPTION  OF  ENVIRONMENTS 

Four  different  environments  were  evaluated.  We  describe  these  according 
to  the  following  general  characteristics.  The  first  two  environments  have  range- 
independent  ocean  bottom  depths,  while  the  third  and  fourth  environments  have 
range-dependent  bottom  depths.  All  environments  have  range-independent 
sound  speed  profiles  with  a  bilinear  shape.  Also,  the  first  three  environments  are 
specified  as  a  single,  angle-independent  radial,  of  maximum  range  equal  to 
15000  m.  The  fourth  environment  is  specified  out  to  40000  m. 

The  first  environment  has  a  constant  600  m  water  depth,  and  an  SSP  with 
a  positive  (downward  refracting)  upper  layer  and  a  negative  (upward  refracting) 
lower  layer.  The  second  environment  has  a  constant  120  m  water  depth,  and  an 
SSP  with  both  layers  having  a  negative  (upward  refracting)  slope.  The  third 
environment  has  a  range-dependent  bottom  depth  of  600  m  at  the  source 
position,  and  ends  with  a  200  m  bottom  depth  at  15000  m.  The  final  environment 
is  intended  to  represent  a  deeper  water  ocean  having  a  3000  m  depth  at  the 
source  position,  and  ending  with  a  500  m  depth  at  40000  m.  These  environments 
are  intended  to  represent  a  variety  of  conditions  such  as  a  very  shallow  water 
littoral,  medium  depth  coastal  shelf,  and  deep  water  transition  region. 
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The  density  of  the  water  in  all  environments  was  set  at  1024  kg/m^.  The 
water  absorption  was  specified  by  the  Fisher-Simmons  model  where  the  salinity 
=  35  ppm,  temperature  =  15.6  °C,  and  pH  =  8.  In  all  cases,  the  bottom  was 
specified  as  a  single  fluid  layer  with  no  shear  properties.  The  bottom  parameters 

were  specified  as  Cb  =  1600  m/s,  =  0.2  s'\  pb  =  1700  kg/m^,  ab  =  0.1  dB/A,. 

The  bottom  roughness  was  produced  with  the  Fox-Hayes  approach,  where  a  =  1 
m,  Uor  =  100  m,  and  p  =  3.5  (see  Eq.  (77)).  In  all  cases,  the  frequency  of 
transmission  was  100  Hz.  Also,  the  source  velocity  =  15  kts,  and  the  receiver 
velocity  =  -5  kts  (indicating  a  receiver  traveling  in  the  opposite  direction  of  the 
source).  In  each  numerical  experiment,  the  source  depth  was  run  at  both  =  50 
m  and  =  3  m  in  order  to  compare  the  impact  of  source  depth. 

A  ray  trace  of  each  environment  was  computed  using  the  method 
described  in  Appendix  B.  The  implementation  supports  a  range-dependent  SSP 
and  bottom  depth,  and  assumes  perfect  surface  and  bottom  reflection.  Appendix 
B  includes  a  derivation  of  the  method  along  with  a  User’s  Guide  supporting  its 
operation. 

The  horizontal  and  vertical  mesh  sizes,  along  with  the  frequency 
resolution  associated  with  each  environment,  are  listed  in  Table  4. 


Table  4.  Mesh  sizes  used  for  each  environment 


Environment 

1 

2 

3 

4 

Nr 

512 

512 

512 

512 

Ar 

29.3  m 

29.3  m 

29.3  m 

78.13  m 

Nz 

1024 

512 

1024 

2048 

Az 

1.17  m 

0.47  m 

1.17  m 

2.93  m 

Nf 

60 

60 

60 

60 

Af 

0.0132  Hz 

0.0132  Hz 

0.0132  Hz 

0.0132  Hz 

It  is  recognized  that  the  above  parameters  provide  a  somewhat  courser 
range  and  (in  some  cases)  depth  resolution  than  the  recommendations  given  in 
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Section  III.E.  In  the  first  three  environments,  at  100  Hz,  this  represents  about  a 
factor  of  two  more  than  the  recommended  mesh  sizes  in  Eqs.  (81)  and  (82). 
However,  our  analysis  in  Section  IV  showed  that  the  degradation  is  relatively 
minor  out  to  moderate  ranges  when  the  mesh  size  is  one  half  of  what  might  be 
considered  optimal.  In  our  fourth  environment,  the  mesh  sizes  are  quite  coarser 
than  recommended.  However,  it  is  felt  that  the  resolution  is  sufficiently  adequate 
to  characterize  the  shape  of  the  TL  response  in  frequency  versus  range  and 
depth,  and  that  the  absolute  magnitude  values  are  of  secondary  interest. 
Limitations  on  memory  and  computational  throughput  dominated  our  rationale  for 
the  mesh  size  selected  in  this  case. 

The  acoustic  propagation  for  each  environment  and  the  two  source  depths 
is  described  by  a  ray  trace  plot,  source  ^-/(wavenumber-frequency)  contour  plot, 
and  a  TL  versus  range  and  depth  plot  representing  the  minimum  TL  across 
frequency.  Also,  the  TL  versus  frequency  response  for  various  ranges  and  a 
fixed  depth  was  plotted  for  each  experiment. 

The  bandwidth  (BW)  parameter  described  by  Eq.  (84)  was  computed  for 
each  environment  and  source  depth  case,  and  at  each  range  and  depth  cell.  The 
aggregate  behavior  is  described  by  a  histogram  of  BW  variables  from  all  cells  in 
the  water  column.  The  density  of  BW  is  plotted  along  a  horizontal  axis  that  spans 
a  zero  to  a  six-a  range.  Also,  the  BW  versus  the  normalized  magnitude  of 
relative  pressure  at  each  cell  is  plotted  on  a  2-D  scatter  plot.  The  normalization 

refers  to  multiplying  the  PE  pressure  field  by  Vr  to  remove  the  effects  of 
cylindrical  spreading  on  the  statistics.  Thus,  the  scatter  plots  represent  the 
dependency  of  BW  on  the  local  interference  structure,  vice  the  absolute 
transmission  loss.  Finally,  the  average  BW  across  depth  for  each  range  cell,  and 
across  range  for  each  depth  cell  is  plotted.  These  plots  show  the  variability  of  a 
marginalized  BW  with  range  and  depth. 
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B.  DESCRIPTION  OF  RESULTS 


This  section  describes  the  results  associated  with  each  environment. 

1.  Environment#! 

The  depth  dependent  acoustical  parameters  (SSP,  density  and 
attenuation)  for  the  first  environment  are  illustrated  in  Fig.  27.  A  set  of  ray  traces 
projected  within  ±10°  vertical  angle  limits,  at  a  source  depth  of  50  m,  is  illustrated 
in  Fig.  28.  At  the  source  depth  of  50  m,  the  rays  show  a  combination  of  refracted 
direct  paths  and  bottom  bounce  propagation. 

Speed  of  Sound  (m/s)  Density  (kg/m^)  Attenuation  (dB/Km) 


Figure  27.  Acoustic  parameters  for  environment  #1 


Ray  Trace  Plot,  11  rays,  start  angles  [-10.0,10.0]  degs,  zs  =  50  m 
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Figure  28.  Ray  trace  for  environment  #1 ,  Zs  =  50  m 
Figure  29  is  a  contour  plot  of  the  distribution  of  source  magnitude  across 
wavenumber  and  frequency.  The  figure  represents  the  distribution  of  the  energy 
caused  by  a  15  knot  source  speed  at  100  FIz,  and  a  0°  vertical  direction.  The 
original  CW  source  is  transmitted  into  the  medium  with  a  min-max  bandwidth  of 
about  0.52  Hz  due  to  source  motion-induced  Doppler  spreading.  The  magnitudes 
of  the  source  function  vary  with  wavenumber  due  to  the  Thomson-Bohun^^  taper 
function  described  previously  in  Eq.  32. 

Transmitted  source  vs  k-f 


100  100.1  100.2  100.3  100.4  100.5  100.6  100.7 

frequency  {Hz) 

Figure  29.  k-f  spectrum  of  source  in  environment  #  1 ,  Zs  =  50  m,  Vs  =  1 5  kts 

Figure  30  illustrates  the  range  and  depth-dependent  transmission  loss 
computed  by  the  M3PE  implementation  for  environment  #1  out  to  15,000  m  with 
the  source  at  50  m.  The  plot  represents  the  minimum  loss  at  each  range-depth 
cell  for  all  frequencies  evaluated.  It  can  be  observed  that  dominant  energy  is 
propagated  at  fairly  shallow  angles.  High  angle  propagation  is  also  apparent  in 
the  interference  pattern,  but  with  higher  transmission  loss.  Figure  31  plots  a  2-D 
slice  of  the  TL  vs  range  at  the  depth  of  50  m. 
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Transmission  loss  (dB).  f  =  100  Hz,  =  50  m,  =  15  Kts,  v^  =  -5  Kts 
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Figure  30.  Minimum  transmission  loss  across  frequency  vs  range  and  depth  for 

environment  #  1 ,  Zs  =  50  m,  Vs  =  15  kts,  Vr  =  -5  kts 


Transmission  Loss  vs  range,  =  50  m,  f  =  1 00  Hz.  =  50  m.  =  1 5  Kts,  =  -5  Kts 


Figure  31 .  Minimum  transmission  loss  across  frequency  vs  range  for  environment 

#  1,  Zs  =  50  m,  Zr  =  50m 
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Figure  32  is  an  example  plot  of  the  transmission  loss  at  all  frequencies, 
evaluated  at  the  discrete  ranges  between  5000  m  and  5300  m,  at  a  depth  of  50 
m.  From  the  range  resolutions  provided  in  Table  4,  this  translates  into  about  10 
separate  response  curves  overlaid  on  the  same  plot.  Interestingly,  it  can  be 
observed  that  most  of  the  energy  exists  above  the  maximum  transmitted 
frequency  indicated  by  Fig.  29.  This  is  because  the  receiver  motion  (-5  knots),  in 
the  opposite  direction  of  the  source  motion,  imparts  an  extra  Doppler  shift  to  the 
TL  distribution  towards  the  right.  Thus,  the  dominant  energy  now  exists  about  a 
frequency  of  100.65  Hz.  It  can  also  be  observed  that  the  main  lobe  has  a  width  of 
about  0.08  Hz  at  the  10  dB  down  points.  Also,  there  is  significant  “structure”  at 
other  frequencies  (namely  a  set  of  local  maxima  around  100.55  Hz)  although  at  a 
much  lower  level  (~40  dB  below  the  peak). 


Transmission  Loss  (dB)  vs  frequency,  =  50  m,  v^  =  15  Kts,  v^  =  -5  Kts 


frequency  (Hz) 

Figure  32.  Transmission  loss  vs  frequency,  environment  #1 ,  Zs  =  50  m,  Zr  =  50  m, 

5000  m  <  r  <  5300  m,  Vs  =  15  kts,  Vr  =  -5  kts 

Figure  33  illustrates  TL  distribution  across  frequency  for  all  ranges 
between  0-15000  m,  at  the  constant  received  depth  of  50  m.  In  this  3-D  plot,  it 
can  be  seen  that  the  TL  response  is  fairly  broadband  at  short  ranges,  but  rapidly 
tapers  to  a  steady  state  narrow  shape  after  about  1-2  km.  This  near-field  broad 
bandwidth  is  probably  caused  by  the  predominance  of  high  angle 
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propagation  that  still  has  significant  energy  close  to  the  source.  This  effect  is 
likely  also  amplified  by  the  Thomson  and  Bohun  high  wavenumber  correction 
factor  incorporated  into  the  source  function.  This  energy  rapidly  attenuates  after 
high  angle  boundary  interactions,  causing  the  shape  of  the  TL  response  to 
concentrate  around  a  single  local  maximum  after  a  few  thousand  meters. 


TL  vs  range  and  frequency,  z  =  50  m,  f  =  100  Hz,  z  ^  =  50  m,  =  15  Kts,  =  -5  Kts 


Figure  33.  Transmission  loss  vs  frequency  and  range,  environment  #1 ,  Zs  =  50  m, 

Zr  =  50  m.  Vs  =  1 5  kts,  Vr  =  -5  kts 


The  aggregate  behavior  of  the  bandwidth  (BW)  estimates  at  each  range 
and  depth  cell,  for  environment  1  with  the  source  at  50  m,  is  depicted  in  Figs.  34 
and  35.  Recall  that  the  bandwidth  parameters  are  computed  on  the  linear  PE 
pressure  magnitudes  with  the  frequency-independent  cylindrical  spreading  term 
removed.  Figure  34a  shows  the  average  BW  at  each  range  across  all  depths  in 
the  water  column.  Figure  34b  shows  the  average  BW  at  each  water  depth  across 
all  range  cells.  It  can  be  seen  that  the  BW  depends  somewhat  on  range,  but  not 
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as  much  on  depth  in  this  environment.  Figure  34a  is  simply  a  parametric 
representation  of  the  frequency  response  characteristics  observed  in  the  3-D  TL 
distribution  shown  in  Fig.  33. 


Marginal  Distribution  of  Bandwidth  vs  Range,  =  50  (m),  v^  =  15  Kts,  =  -5  Kts 


MarginAlDistntMJtion  of  Bandwidth  vs  Depth,  7  -50(m).v  -ISKts.v^  5Kts 
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Figure  34.  Marginalized  bandwidth  vs  range  (a),  and  depth  (b) ,  environment  #1 , 

Zs  =  50  m,  Vs  =  15  kts,  Vr  =  -5kts 

Figure  35a  shows  an  overall  histogram  of  bandwidths  for  all  ranges  and 
depths  (above  the  bottom).  From  this  plot,  a  sample  mean  BW  of  about  0.06  FIz, 
and  a  standard  deviation  of  about  0.05  Hz  has  been  computed.  Finally,  the 
dependency  of  BW  on  the  local  peak  pressure  magnitudes  across  frequency  is 
provided  in  Fig.  35b.  At  first  glance,  this  plot  might  seem  to  contradict  the 
behavior  shown  in  Figs.  33  and  34a,  where  shorter  ranges  (and  lower  TL) 
provide  for  larger  bandwidths.  That  is,  the  overall  trend  in  Fig.  35b  is  for  cells  with 
high  pressure  magnitudes  to  have  smaller  bandwidths  than  cells  with  lower 
magnitudes.  However,  a  second  look  at  this  figure  indicates  that  there  is  a 
segment  of  samples  with  high  bandwidth  at  the  larger  pressure  magnitudes. 
These  bandwidths  are  likely  from  the  closer  ranges.  Thus  it  is  concluded  that  two 
contradictory  influences  are  in  play.  We  can  summarize  this  by  stating  that  the 
bandwidth  values  are  generally  proportional  to  absolute  TL  level,  except  at  short 
ranges  where  high  bandwidths  are  caused  by  high  angle  propagation,  even  when 
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the  overall  TL  is  low. 


PDF  of  baixJwidth  for  all  ranges  and  depths  in  water  column,  Zg  =  50  (m),  =  1 5  Kts,  v^  =  -5  Kts 


(a)  (b) 

Figure  35.  Aggregate  histogram  (a)  and  scatter  plot  vs  normalized  linear  TL  (b), 

environment  #1 ,  Zs  =  50  m,  Vs  =  1 5  kts,  Vr  =  -5kts 

The  next  set  of  Figs.  36-42  repeat  the  results  presented  in  Figs.  28-35  for 
the  same  environment  (#1),  but  with  the  source  positioned  at  3  m  instead  of  50 
m.  Here,  the  ray  trace  (Fig.  36)  and  minimum  TL  plots  (Figs.  37  and  38)  show  a 
TL  response  that  is  much  more  dominated  by  high  angle  propagation  and 
boundary  interactions.  The  shape  of  the  frequency  response  shown  in  Figs.  39 
and  40  is  similar,  but  the  dominant  frequency  component  is  a  bit  wider.  Also,  a 
secondary  propagation  mode  at  about  100.37  Hz  is  more  dominant  for  the 
shallow  source.  Figures  41  and  42  show  a  very  similar  behavior  as  for  the  deeper 
source,  except  that  all  plots  are  shifted  towards  higher  bandwidths  by  about  0.02 
Hz.  Thus,  the  shallower  source  demonstrates  a  marginally  higher  aggregate 
bandwidth  than  the  deeper  source  in  this  environment.  An  interesting  difference 
between  the  two  sets  of  results  is  that  the  average  minimum  TL  at  the  outer 
ranges  is  about  5-10  dB  greater  for  the  shallow  source  than  the  deeper  one.  As 
expected,  the  higher  angle  propagation  also  results  in  more  TL. 
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Ray  Trace  Plot,  1 1  rays,  start  angles  [-10.0,10.0]  degs,  zs  =  3  m 


10000 


range  (m) 

Figure  36.  Ray  trace  for  environment  #1 ,  Zs  =  3  m 
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Figure  37.  Minimum  transmission  loss  across  frequency  vs  range  and  depth  for 

environment  #1,Zs  =  3m,Vs  =  15  kts,  Vr  =  -5  kts 
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Transmission  Loss  vs  range,  z  ,=  50  m.f=  100  Hz,  z  =3m,  v  =  15  Kts,  v  =-5  Kts 
^  plot  s  s  r 


Figure  38.  Minimum  transmission  loss  across  frequency  vs  range  for  environment 

#  1 ,  Zs  =  3  m,  Zr  =  50m,  Vs  =  1 5  kts,  Vr  =  -5  kts 


Transmission  Loss  (dB)  vs  frequency,  =  3  m,  =  15  Kts,  v^  =  -5  Kts 


frequency  (Hz) 

Figure  39.  Transmission  loss  vs  frequency,  environment  #1 ,  Zs  =  3  m,  Zr  =  50  m, 

5000  m  <  r  <  5300  m,  Vs  =  15  kts,  Vr  =  -5  kts 
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range  (m) 


TL  vs  range  and  frequency,!  p|„,  =  50  m,f  =  100  Hz. z  ^  =  3  m,  =  15  Kts,  v^  =  -5  Kts 


Figure  40.  Transmission  loss  vs  frequency  and  range,  environment  #1 ,  Zs  =  3  m, 

Zr  =  50  m,  Vs  =  1 5  kts,  Vr  =  -5  kts 


Marginal  Distribution  of  Bandwidth  vs  Range,  z  =3(m).v  =  15  Kts.  v  = -5  Kts 
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Figure  41 .  Marginalized  bandwidth  vs  range  (a),  and  depth  (b) ,  environment  #1 , 
Zs  =  3  m.  Vs  =  1 5  kts,  Vr  =  -  5kts 
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PDF  of  bandwidth  for  all  ranges  and  depths  in  water  column.  Zg  =  3  (m),  =  1 5  Kts, 


Distribution  of  Bandwidth  vs  Pressure  Magnitude,!  ^  =  3  (m),  =  15  Kts,  v^  =  -5  Kts 


0.025 

0.02 

■o 


|>  0.015 


bandwidth  (Hz) 


(a)  (b) 

Figure  42.  Aggregate  histogram  (a)  and  scatter  plot  vs  normalized  linear  TL  (b), 

environment  #1 ,  Zs  =  3  m,  Vs  =  15  kts,  Vr  =  -5kts 


2.  Environment  #2 

The  depth-dependent  parameters  associated  with  environment  #2  are 
illustrated  in  Fig.  43.  The  ray  trace  in  Fig.  44  indicates  that  this  environment  is  a 
very  shallow  water  waveguide  with  upward  refracting  propagation.  The  source 
spectrum  versus  wavenumber  is  shown  in  Fig.  45.  The  transmission  loss  results 
versus  range  and  depth  for  a  source  at  50  m  are  shown  in  Figs.  46-47.  It  can  be 
observed  that  relatively  low  TL  values  occur  in  the  surface  duct  even  at  long 
ranges  when  the  source  is  deep.  This  is  confirmed  by  the  ray  trace  which  shows 
that  propagation  at  narrow  angles  from  a  source  at  50  m  can  occur  with  large 
skip  distances  between  surface  interactions.  The  TL  versus  frequency  and  range 
is  shown  in  Figs.  48  and  49.  As  a  side  note,  Fig.  46  also  allows  one  to  observe 
an  example  of  the  rough  bottom  generated  for  all  environments.  In  this  plot  (vice 
the  other  environments),  the  ocean  is  sufficiently  shallow  such  that  the  scale 
allows  sufficient  bottom  structure  to  be  observable. 

The  bandwidth  parameters  are  illustrated  in  Figs.  50  and  51.  The  mean 
bandwidth  is  about  0.06  Hz  and  the  standard  deviation  was  0.03  Hz.  Figures  49 
and  50a  illustrate  that  the  transition  between  wide  bandwidths  at  short  range,  and 
the  steady  state  narrow  bandwidths,  which  occur  much  more  quickly  in 
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environment  #2  than  environment  #1 .  However,  the  dependence  of  bandwidth  on 
receiver  depth  is  more  significant  in  this  environment  than  the  previous  one. 
Specifically,  the  bandwidths  are  wider  when  the  receiver  is  near  the  boundaries, 
vice  towards  the  middle  of  the  water  column.  This  is  corroborated  by  the  fact  that 
there  is  generally  more  TL  near  the  boundaries  than  in  the  middle  of  the 
waveguide  in  this  environment.  Figure  51b  also  indicates  that  the  bandwidths  are 
generally  wider  when  the  pressure  magnitude  is  lower  (i.e.,  higher  TL). 

Speed  of  Sound  (m/s)  Density  (kg/m^)  Attenuation  (dB/Km) 
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Figure  43.  Acoustic  parameters  for  environment  #2 


Figure  44.  Ray  trace  for  environment  #2,  Zs  =  50  m 
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Transmitted  source  vs  k-f 


Figure  45.  k-f  spectrum  of  source  in  environment  #  2,  Zs  =  50  m,  Vs  =  15  kts 
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Figure  46.  Minimum  transmission  loss  across  frequency  vs  range  and  depth  for 

environment  #  2,  Zs  =  50  m,  Vs  =  15  kts,  Vr  =  -5  kts 
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Transmission  Loss  vs  range,  =  50  m,  f  =  100  Hz,  =  50  m.  =  15  Kts,  =  -5  Kts 


Figure  47.  Minimum  transmission  loss  across  frequency  vs  range  for  environment 

#  2,  Zs  =  50  m,  Zr  =  50m,  Vs  =  15  kts,  Vr  =  -5  kts 


Transmission  Loss  (dB)  vs  frequency,  =  50  m,  =  15  Kts,  v^  =  -5  Kts 


Figure  48.  Transmission  loss  vs  frequency,  environment  #2,  Zs  =  50  m,  Zr  =  50  m, 

5000  m  <  r  <  5300  m,  Vs  =  15  kts,  Vr  =  -5  kts 
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Figure  49.  Transmission  loss  vs  frequency  and  range,  environment  #2,  Zs  =  50  m, 

Zr  =  50  m,  Vs  =  1 5  kts,  Vr  =  -5  kts 
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Marginal  Distribution  of  Bandwidth  vs  Depth,  z^  =  50  (m),  =  1 5  Kts,  v^  =  -5  Kts 
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Figure  50.  Marginalized  bandwidth  vs  range  (a),  and  depth  (b),  environment  #2,  Zs 

=  50  m.  Vs  =  15  kts,  Vr  =  -5kts 
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PDF  of  bandwidth  for  all  ranges  and  depths  in  water  column,  =  50  (m),  =  15  Kts,  = 


Distribution  of  Bandwidth  vs  Pressure  Magnitude,  z  ^  =  50  (m),  =  15  Kts,  =  -5  Kts 


0  0.02  0,04  0.06  0.08  0,1  0.12  0.14  0.16  0.18 

bandwidth  (Hz) 


(a)  (b) 

Figure  51 .  Aggregate  histogram  (a)  and  scatter  plot  vs  normalized  linear  TL  (b), 

environment  #2,  Zs  =  50  m,  Vg  =  15  kts,  Vr  =  -5kts 

Figures  52-58  repeat  the  numerical  experiments  associated  with 
environment  #2  for  a  source  at  3  m.  Again  we  observe  (in  Figs.  53  and  54)  that  a 
source  at  the  shallow  depth  has  generally  more  TL  (at  comparable  ranges)  than 
a  deeper  source.  Also,  Figs.  55  and  56  indicate  that  more  relative  energy  is 
distributed  across  frequency  when  the  source  depth  is  shallow.  Interestingly,  the 
secondary  contributions  at  the  lower  frequencies,  for  the  3  m  source  depth,  are 
only  about  12  dB  down  from  the  main  peak,  while  they  are  almost  40  dB  down  for 
the  50  m  source  depth  case.  The  average  BW  estimates  plotted  in  Fig.  57(a,b) 
with  a  source  at  3  m  are  generally  about  0.05  Hz  greater  than  the  BW  values 
from  a  source  at  50  m,  for  environment  #2.  This  could  again  be  due  to  a  greater 
tendency  toward  high  angle  propagation,  as  indicated  by  the  TL  versus  range 
and  depth  plot  in  Fig.  53.  The  results  for  environments  #1  and  #2  demonstrate 
how  two  very  different  sets  of  acoustic  conditions  (i.e.,  one  with  downward 
refraction  and  the  other  with  upward  refraction)  can  still  have  similar  effects  on 
the  characteristics  of  received  signal  bandwidth. 

Also,  the  relationship  between  bandwidth  and  normalized  pressure, 
plotted  in  Fig.  58b,  is  less  dispersed  for  the  shallow  source  than  the  deeper  one. 
This  is  likely  due  to  the  fact  that  the  range  normalized  transmission  loss  is  less 
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variable  (i.e.,  more  uniform  across  depth)  when  the  source  is  shallow  than  for  the 
deeper  source  in  this  particular  environment.  This  is  also  supported  by  Fig.  57b 
which  indicates  very  consistent  bandwidth  parameters  across  depth,  in  contrast 
to  the  plot  in  Fig.  50b  for  the  deep  source. 


Figure  52.  Ray  trace  for  environment  #2,  Zs  =  3  m 
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Figure  53.  Minimum  transmission  loss  across  frequency  vs  range  and  depth  for 

environment  #  2,  Zs  =  3  m,  Vg  =  15  kts,  Vr  =  -5kts 


Figure  54.  Minimum  transmission  loss  across  frequency  vs  range  for  environment 

#  2,  Zg  =  3  m,  Zr  =  50m,  Vg  =  15  kts,  Vr  =  -5kts 


Transmission  Loss  (dB)  vs  frequency,  =  3  m,  =  1 5  Kts,  =  -5  Kts 


Figure  55.  Transmission  loss  vs  frequency,  environment  #2,  Zg  =  3  m,  Zr  =  50  m, 

5000  m  <  r  <  5300  m,  Vg  =  15  kts,  Vr  =  -5kts 
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range  (m) 


TL  vs  range  and  frequency,  z  =  50  m.f  =  100  Hz.i  ^  =  3  m,  =  15  Kts,  v^  =  -5  Kts 


Figure  56.  Transmission  loss  vs  frequency  and  range,  environment  #2,  Zg  =  3  m, 

Zr  =  50  m,  Vg  =  15  kts,  Vr  =  -5kts 


Marginal  Distribution  of  Bandwidth  vs  Range.  Zg  =  3  (m),  =  15  Kts.  v^  = 


Marginal  Distribution  of  Bandwidth  vs  Depth,  =  3  (m),  =  15  Kts,  v^  =  -5  Kts 


Figure  57. 


(a)  (b) 

Marginalized  bandwidth  vs  range  (a),  and  depth  (b) , 
Zg  =  3  m,  Vg  =  15  kts,  Vr  =  -5kts 


environment  #2, 
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Distribution  of  Bandwidth  vs  Pressure  Magnitude,!  ^  =  3  (m),Vg  =  15  Kts,v^  =  -5  Kts 


PDF  of  bandwidth  for  all  ranges  and  depths  in  water  column,  =  3  (m),  =  1 5  Kts.  = 
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Figure  58. 
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Aggregate  histogram  (a)  and  scatter  plot  vs  normalized  linear  TL  (b) 
environment  #2,  Zs  =  3  m,  Vs  =  15  kts,  Vr  =  -5kts 


3.  Environment  #3 

The  next  two  environments  demonstrate  the  characteristics  of  received 
bandwidth  when  the  ocean  depth  is  range-dependent.  The  depth-dependent 
acoustic  parameters  for  environment  #3  are  plotted  in  Fig.  59.  Figures  60-67 
contain  results  when  the  source  depth  is  50  m.  Figures  68-74  provide  the  results 
when  the  source  depth  is  3  m.  The  ray  trace  in  Fig.  60  shows  a  mini  convergence 
zone  (with  some  bottom  bounce  at  steeper  angles)  in  the  first  8-9  km,  followed  by 
multiple  bottom  bounces  in  the  last  5  km  when  the  ocean  depth  becomes 
shallow.  Figure  68  demonstrates  bottom  bounce  propagation  at  smaller  launch 
angles  when  the  source  is  shallow.  The  TL  versus  range  and  depth  plot  in  Fig.  62 
again  indicates  fairly  shallow  angle  propagation  when  the  source  is  deep.  This 
contrasts  with  the  much  steeper  propagation  illustrated  in  Fig.  69  in  the  first 
10000  m. 

The  frequency  versus  range  TL  response  given  in  Figs.  64-65  (deep 
source)  and  71-72  (shallow  source)  possess  a  wider  main  component  and  more 
consistent  broadband  pedestal  when  the  source  is  shallow  than  when  it  is  deep. 
In  Fig.  67a,  the  mean  BW  for  the  deep  source  has  been  computed  to  be  about 
0.06  Hz.  When  the  source  is  shallow.  Fig.  74a  indicates  a  slightly  larger  mean 
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BW  of  about  0.09  Hz,  consistent  with  the  behavior  observed  in  the  first  two 
environments.  The  scatter  plots  in  Figs.  67b  and  74b  indicate  a  unimodal 
dependence  upon  local  magnitude  when  the  source  is  deep,  but  a  somewhat 
bimodal  structure  when  the  source  is  shallow.  This  will  be  discussed  in  more 
detail  below. 

A  unique  characteristic  of  this  range-dependent  environment  is  related  to 
a  fairly  substantial  dependence  of  the  BW  parameter  on  range  and  depth.  This 
can  be  observed  in  Figs.  66  and  73,  regardless  of  the  source  depth.  Specifically, 
the  BW  undergoes  a  significant  decrease  after  a  range  of  about  10000  m,  and 
above  a  depth  of  200  m.  This  corresponds  to  propagation  above  the  shallow 
plateau  during  the  last  5000  m  of  range.  The  TL  responses  in  Figs.  62  and  69 
indicate  that  acoustic  propagation  above  this  shelf  occurs  with  fairly  low  grazing 
angles,  even  though  the  ray  trace  in  Fig.  60  might  suggest  otherwise.  It  is 
possible  that  propagation  angles  increased  past  the  critical  angle  once  the  wave- 
front  passed  over  the  shallow  shelf,  and  thus  rapidly  attenuated.  Thus,  in  this 
environment,  the  bottom  profile  actually  appears  to  filter  high  angle  propagation 
at  the  long  ranges  causing  the  bandwidth  to  decrease  during  this  latter  segment 
of  range.  The  differences  noted  in  the  mean  BW  (across  all  range  and  depth 
cells),  in  Figs.  67  and  74,  are  likely  due  to  the  preponderance  of  high  angle 
propagation  in  the  first  10000  m  that  occurs  when  the  source  is  shallow. 

The  marked  difference  in  propagation  angles  during  the  first  10000  m  and 
the  latter  5000  m  when  the  source  is  shallow  is  probably  the  cause  of  the 
bimodality  in  the  histogram  and  scatter  plots  noted  above.  When  the  source  is 
deep,  the  propagation  angles  are  more  consistent  over  range,  resulting  in  a  more 
unimodal  density  function.  This  numerical  experiment  has  allowed  us  to  observe 
several  interesting  impacts  on  the  characteristics  of  the  bandwidth 
parameter,when  the  environment  has  more  complicated  range-dependent 
attributes. 
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Figure  59.  Acoustic  parameters  for  environment  #3 


Ray  Trace  Plot,  1 1  rays,  start  angles  (-10  0,10  0]  degs,  zs  =  50  m 
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Figure  60.  Ray  trace  for  environment  #3,  Zs  =  50  m 
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Transmitted  source  vs  k-f 


frequency  (Hz) 


Figure  61 .  k-f  spectrum  of  source  in  environment  #  3,  Zs  =  50  m,  Vs  =  1 5  kts 


range  (m) 

Figure  62.  Minimum  transmission  loss  across  frequency  vs  range  and  depth  for 

environment  #  3,  Zs  =  50  m,  Vs  =  15  kts,  Vr  =  -5  kts 
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Transmission  Loss  vs  range,  =  50  m,  f  =  100  Hz,  =  50  m,  v^  =  15  Kts,  =  -5  Kts 


Figure  63.  Minimum  transmission  loss  across  frequency  vs  range  for  environment 

#  3,  Zs  =  50  m,  Zr  =  50m,  Vs  =  15  kts,  Vr  =  -5kts 


Transmission  Loss  (dB)  vs  frequency,  z^  =  50  m.  v^  =  15  Kts,  v^  =  *5  Kts 


Figure  64.  Transmission  loss  vs  frequency,  environment  #3,  Zs  =  50  m,  Zr  =  50  m, 

5000  m  <  r  <  5300  m,  Vs  =  15  kts,  Vr  =  -5  kts 
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TL  vs  range  and  frequency,  z  =  50  m.f  =  100  Hz,  z  ^  =  50  m,  =  15  Kts,  v^  =  -5  Kts 
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Figure  65.  Transmission  loss  vs  frequency  and  range,  environment  #3,  Zs  =  50  m, 

Zr  =  50  m,  Vs  =  1 5  kts,  Vr  =  -5  kts 


(a)  (b) 

Figure  66.  Marginalized  bandwidth  vs  range  (a),  and  depth  (b),  environment  #3,  Zs 

=  50  m.  Vs  =  15  kts,  Vr  =  -5kts 
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PDF  of  bandwidth  for  all  ranges  and  depths  in  water  column,  =  50  (m).  =  1 5  Kts, 


Distribution  of  Bandwidth  vs  Pressure  Magnitude,!  ^  =  50  (m),  =  15  Kts,  =  -5  Kts 
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Figure  67.  Aggregate  histogram  (a)  and  scatter  plot  vs  normalized  linear  TL  (b), 

environment  #3,  Zs  =  50  m,  Vs  =  15  kts,  Vr  =  -5kts 


RayTrace  Plot,  11  rays,  start  angles  [-10.0,10.0]  degs,zs  =  3  m 
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Figure  68.  Ray  trace  for  environment  #3,  Zs  =  3  m 
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Figure  69.  Minimum  transmission  loss  across  frequency  vs  range  and  depth  for 

environment  #  3,  Zs  =  3  m,  Vs  =  15  kts,  Vr  =  -5kts 


Transmission  Loss  vs  range,  z  ,=  50  m,  f=  100  Hz,  z  =3m,v  =  15  Kts,  v  = -5  Kts 
^  plot  ’  s  ’  s  ’  r 


Figure  70.  Minimum  transmission  loss  across  frequency  vs  range  for  environment 

#  3,  Zs  =  3  m,  Zr  =  50m,  Vs  =  15  kts,  Vr  =  -5kts 


80 


Transmission  Loss  (dB)  vs  frequency,  =  3  m,  =  1 5  Kts,  v^  =  -5  Kts 


Figure  71 .  Transmission  loss  vs  frequency,  environment  #3,  Zs  =  3  m,  Zr  =  50  m, 

5000  m  <  r  <  5300  m,  Vs  =  15  kts,  Vr  =  -5  kts 


TL  vs  range  and  frequency,  z  =  50  m.f  =  100  Hz,  z  ^  =  3  m,  =  15  Kts,  v^  =  -5  Kts 


Figure  72.  Transmission  loss  vs  frequency  and  range,  environment  #3,  Zg  =  3  m, 

Zr  =  50  m.  Vs  =  1 5  kts,  Vr  =  -5  kts 
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Marginal  Distribution  of  Bandwidth  vs  Range,  Zg  =  3  (m),  v^  =  1 5  Kts,  =  -5  h 


Marginal  Distribution  of  Bandwidth  vs  Depth,  Zg  =  3  (m),  v^  =  15  Kts,  v^  =  -5  Kts 
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Figure  73.  Marginalized  bandwidth  vs  range  (a),  and  depth  (b),  environment  #3,  Zg 

=  3  m,  Vg  =  15  kts,  Vr  =  -5kts 
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Figure  74.  Aggregate  histogram  (a)  and  scatter  plot  vs  normalized  linear  TL  (b), 

environment  #3,  Zg  =  3  m,  Vg  =  15  kts,  Vr  =  -5kts 


4.  Environment  #4 

The  final  environment  is  designed  to  demonstrate  the  characteristics  of 
received  bandwidth  in  a  fairly  deep,  but  range-dependent  ocean,  such  as  may  be 
encountered  near  shelf  break  regions.  The  depth  dependent  acoustic  parameters 
for  this  environment  are  illustrated  in  Fig.  75.  A  ray  trace  from  a  source  depth  of 
50  m,  depicted  in  Fig.  76,  shows  convergence  zone  (CZ)  propagation  in  the  first 
30000  m,  followed  by  multiple  bottom  bounces  above  the  shallow  shelf.  When 
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the  source  is  shallow,  Fig.  84  indicates  that  the  CZ  is  augmented  with  a  deep 
bottom  bounce  path.  The  minimum  TL  response,  given  in  Figs.  77  (for  the  deep 
source),  and  85  (for  the  shallow  source),  are  consistent  with  the  ray  traces. 
Similar  to  the  other  environments,  propagation  with  higher  vertical  angles  can  be 
observed  when  the  source  is  shallow  vice  deep.  Figures  78  and  86  illustrate  a 
fairly  deep  shadow  zone  at  the  shallow  depths  between  a  range  of  about  5000 
and  18000  m  for  both  the  deep  and  shallow  source  cases. 

The  transmission  loss  versus  frequency  and  range  is  presented  in  Figs.  81 
(for  the  deep  source)  and  88  (for  the  shallow  source).  In  this  deep  ocean 
environment,  these  plots  are  very  similar.  They  show  significant  structure  at  the 
short  ranges,  followed  by  a  pronounced  dip  in  magnitude  at  the  shadow  zone; 
and  then  a  resurgence  of  narrowband  structure  at  the  longer  ranges.  The  level 
versus  frequency  plots  in  Figs.  80  and  87  (for  the  deep  and  shallow  sources, 
respectively)  illustrate  a  flat  TL  response  within  the  shadow  zone. 

The  bandwidth  estimates  versus  range  and  depth  are  plotted  in  Figs.  82 
for  the  deep  source.  A  similar  set  of  plots  is  contained  in  Fig.  89  for  the  shallow 
source.  Comparison  of  Figs.  83a  and  90a  show  that  the  mean  bandwidth  is 
approximately  0.07  Hz  when  the  source  is  shallow,  and  0.05  Hz  when  the  source 
is  deep.  Thus,  the  shallow  source  is  again  associated  with  an  increase  in 
bandwidth,  although  the  increase  is  modest  in  this  deep  water  environment.  The 
scatter  plots  in  Figs.  83b  and  90b  also  illustrate  significant  differences  in  the 
dependency  of  BW  on  normalized  pressure  magnitude.  Although  it  is  difficult  to 
attribute  a  physical  meaning  to  this  characteristic,  it  is  possible  that  the 
differences  could  be  exploited  by  an  automatic  classification  scheme. 

An  interesting  feature  of  the  marginalized  plots  (Figs.  82  and  89)  is  that 
the  BW  tends  to  decrease  at  the  longer  ranges,  while  increasing  at  shallower 
depths.  This  behavior  is  counter-intuitive  because  the  bottom  profile  requires  that 
the  shallower  depths  occur  at  long  ranges.  The  reason  for  this  can  be  understood 
by  observing  a  color  coded  presentation  of  bandwidth  for  all  ranges  and  depths, 
produced  for  the  case  when  the  source  is  at  3  m.  This  plot,  in  Fig.  91,  shows  that 
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a  concentration  of  high  bandwidth  occurs  in  the  entire  shadow  zone,  consistent 
with  the  plots  shown  in  Figs.  80  and  87.  Thus,  average  bandwidth  over  all  depths 
versus  range  will  tend  to  be  higher  at  the  shorter  ranges.  However,  when 
averaging  over  range,  the  bandwidth  will  be  higher  at  the  shallower  depths.  This 
also  shows  that  even  though  the  model  predicts  significant  bandwidths  in 
particular  parts  of  the  ocean,  sometimes  the  TL  is  so  high  that  signals  with  this 
bandwidth  will  only  be  observable  when  the  levels  of  the  radiating  source  are 
extremely  high. 


Speed  of  Sound  (m/s)  Density  (kg/m^)  Attenuation  (dB/Km) 


Figure  75.  Acoustic  parameters  for  environment  #4 
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Ray  Trace  Plot,  11  rays,  start  angles  [-10.0,10.0]  degs,  zs  =  50m 
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Figure  76.  Ray  trace  for  environment  #4,  Zs  =  50  m 


Transmitted  source  vs  k-f 
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Figure  77.  k-f  spectrum  of  source  in  environment  #  4,  Zs  =  50  m,  Vs  =  15  kts 
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Transmission  loss  (dB),  f  =  100  Hz.  z  =  50  m,  v  =15  Kts.  v  =  -5  Kts 
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Figure  78.  Minimum  transmission  loss  across  frequency  vs  range  and  depth  for 

environment  #  4,  Zs  =  50  m,  Vs  =  15  kts,  Vr  =  -5  kts 


Figure  79.  Minimum  transmission  loss  across  frequency  vs  range  for  environment 

#  4,  Zs  =  50  m,  Zr  =  50m,  Vs  =  15  kts,  Vr  =  -5kts 
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Transmission  Loss  (dB)  vs  frequency,  =  50  m,  =  15  Kts,  v^  =  -5  Kts 


Figure  80.  Transmission  loss  vs  frequency,  environment  #4,  Zs  =  50  m,  Zr  =  50  m, 

5000  m  <  r  <  5300  m,  Vg  =  15  kts,  Vr  =  -5  kts 


TL  vs  range  and  frequency,  z  =  50  m.  f  =  100  Hz.  z  ^  =  50  m. 
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Figure  81 .  Transmission  loss  vs  frequency  and  range,  environment  #4,  Zs  =  50  m, 

Zr  =  50  m.  Vs  =  1 5  kts,  Vr  =  -5  kts 
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Marginal  Distribution  of  Bandwidth  vs  Depth,  =  50  (m),  =  1 5  Kts,  v^.  =  -5  Kts 
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Figure  82.  Marginalized  bandwidth  vs  range  (a),  and  depth  (b),  environment  #4,  Zg 

=  50  m,  Vg  =  15  kts,  Vr  =  -5kts 


PDF  of  bandwidth  for  all  ranges  and  depths  in  water  column,  =  50  (m),  =  1 5  Kts,  =  -5 


Distribution  of  Bandwidth  vs  Pressure  Magnitude,!  ^  =  50  (m),  =  15  Kts,  v^  =  -5  Kts 
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Figure  83.  Aggregate  histogram  (a)  and  scatter  plot  vs  normalized  linear  TL  (b), 

environment  #4,  Zg  =  50  m,  Vg  =  15  kts,  Vr  =  -5kts 
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Figure  85.  Minimum  transmission  loss  across  frequency  vs  range  and  depth  for 

environment  #  4,  Zs  =  3  m,  Vs  =  15  kts,  Vr  =  -5kts 
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Transmission Lossvs range. z  ,=  50  m.f=  100 Hz, z  =3m,v  =  15  Kts, v  =-5 Kts 
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Figure  86.  Minimum  transmission  loss  across  frequency  vs  range  for  environment 

#  4,  Zs  =  3  m,  Zr  =  50m,  Vs  =  15  kts,  Vr  =  -5kts 


Transmission  Loss  (dB)  vs  frequency,  z^  =  3  m,  =  15  Kts,  v^  =  -5  Kts 


frequency  (Hz) 

Figure  87.  Transmission  loss  vs  frequency,  environment  #4,  Zs  =  3  m,  Zr  =  50  m, 

5000  m  <  r  <  5300  m,  Vs  =  15  kts,  Vr  =  -5  kts 
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Figure  88.  Transmission  loss  vs  frequency  and  range,  environment  #4,  Zs  =  3  m, 

Zr  =  50  m,  Vs  =  1 5  kts,  Vr  =  -5  kts 
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Marginal  Distribution  of  Bandwidth  vs  Depth,  Zg  =  3  (m),  v^  =  1 5  Kts,  v^  =  -5  Kts 
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Figure  89.  Marginalized  bandwidth  vs  range  (a),  and  depth  (b),  environment  #4,  Zs 

=  3  m,  Vs  =  15  kts,  Vr  =  -5kts 


91 


PDF  of  bandwidth  for  all  ranges  and  depths  in  water  column,  Zg  =  3  (m),  =  15  Kts,  = 


-5Kt 


Distribution  of  Bandwidth  vs  Pressure  Magnitude,!  ^  =  3  (m),  v^  =  15  Kts,  v^  =  -5  Kts 
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Figure  90.  Aggregate  histogram  (a)  and  scatter  plot  vs  normalized  linear  TL  (b), 

environment  #4,  Zs  =  3  m,  Vs  =  15  kts,  Vr  =  -5kts 


Bandwidth  vs  range  and  depth 


Figure  91 .  Color  coded  plot  of  bandwidth  vs  range  and  depth  for  environment  #4, 

Zs  =  3  m,  Vs  =  15  kts,  Vr  =  -5kts 
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VII.  CONCLUSIONS 


We  have  developed  an  implementation  of  a  PE  acoustic  propagation 
model,  including  a  set  of  metrics  and  statistical  tools,  to  help  us  understand  the 
characteristics  of  the  bandwidth  of  a  CW  signal  when  transmitted  from  a  moving 
source  and  received  by  a  moving  receiver.  Plots  of  transmission  loss  versus 
frequency  and  range  produced  by  the  model  are  an  instructive  way  to  observe 
these  affects  across  the  entire  environment.  The  automatic  method  of  estimating 
bandwidth  and  its  marginals  across  range  and  depth  demonstrated  consistent 
results  that  were  always  explainable  after  manual  analysis  of  the  raw  data. 

Executing  the  model  on  several  diverse  environments  (defined  by  various 
acoustical  parameters)  with  both  deep  and  shallow  source  depths  allowed 
observation  of  several  trends.  First,  it  was  generally  the  case  that  locations  within 
the  medium  corresponding  to  higher  transmission  loss  were  also  correlated  with 
higher  bandwidths.  This  appears  to  be  related  to  the  premise  that  as  the  TL 
increases  for  the  dominant  paths,  the  ordinarily  lower  amplitude  secondary  paths 
(which  tend  to  travel  with  higher  angles  of  propagation)  become  more  significant 
in  a  relative  sense.  A  review  of  the  differences  in  the  main  lobe  of  TL  versus 
frequency  plotted  in  Figs.  48  and  55  demonstrates  this. 

Next,  we  observed  that  it  is  possible  for  the  shape  of  the  range-dependent 
ocean  depth  profile  to  affect  the  tendency  towards  shallow  or  high  angle 
propagation,  and  thereby  influence  the  observed  bandwidth  values.  For  example, 
an  ocean  that  transitions  from  deep  to  shallow  water  depths  can  cause  higher 
angle  modes,  propagating  towards  the  shallow  end,  to  exceed  the  critical  angle 
and  attenuate  rapidly.  In  this  case,  the  bandwidth  marginals  versus  range  or 
depth  will  significantly  decrease  at  the  transition  point.  This  situation  can  also  be 
accompanied  by  a  bimodality  in  the  bandwidth  versus  pressure  magnitude 
scatter  plots  (where  pressure  is  inversely  related  to  TL).  Here,  two  competing 
influences  are  at  work.  Longer  ranges  tend  to  have  higher  TL  and  therefore 
higher  bandwidths.  However,  in  ocean  environments  similar  to  #2  and  #4, 
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defined  in  Section  VI,  the  tendency  towards  high  angle  propagation  can  be 
inversely  proportional  to  the  range.  Thus  we  require  the  numerical  model  to  sort 
out  which  type  of  influence  will  dominate. 

Finally,  we  observed  in  all  the  environments  evaluated  that  a  source  near 
the  surface  will  have  a  higher  overall  bandwidth  than  a  source  located  many 
meters  below  the  surface.  Before  attaching  significance  to  this,  it  is  important  to 
note  that  we  have  not  yet  studied  this  influence  in  its  entirety,  since  only  two 
different  depths  were  evaluated  in  our  examples.  For  example,  it  may  be  possible 
that  a  source  near  the  bottom  produces  bandwidths  that  are  similar  to  those 
produced  by  sources  near  the  surface.  We  also  note  that  only  two  types  of 
(bilinear)  SSPs  were  evaluated,  and  that  the  shape  of  the  SSP  will  likely  have  a 
significant  influence  on  the  bandwidth.  In  the  examples  studied  in  this  thesis,  the 
depth  of  the  source  had  a  considerable  impact  on  the  likelihood  of  high  angle 
propagation,  due  to  the  competing  influences  of  higher  sound  speeds  near  the 
boundaries.  It  remains  to  be  seen  how  robust  this  affect  will  be  when  the  model  is 
exercised  with  actual  SSPs  and  other  acoustical  parameters  corresponding  to 
real  oceans.  It  is  also  unclear  at  this  time  how  rough  sea  surface  scattering  may 
affect  these  results. 

However,  it  is  apparent  from  this  thesis  that  the  source  depth  can  have  a 
major  influence  on  the  dominant  angles  of  propagation  and,  as  a  consequence, 
bandwidth  values.  It  is  then  reasonable  to  conclude  that  information  exists  in  the 
measured  bandwidth  of  a  received  source  that  might  be  helpful  in  determining 
the  likely  depth  of  the  transmitting  source.  Obviously,  BW  effects  could  be 
significantly  useful  in  the  design  of  automatic  classification  schemes  attempting 
to  discriminate  between  sources  on  or  below  the  ocean  surface. 
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APPENDIX  A.  USER’S  GUIDE  FOR  THE  M3PE  MODEL 


This  appendix  provides  a  guide  for  executing  the  M3PE  model  with  user 
defined  inputs  that  define  operating  modes  and  environmental  acoustic 
parameters.  The  first  sub-section  describes  the  various  MATLAB  files  for 
execution,  control,  and  environment  definition.  A  flow  diagram  indicating  the 
functions  performed  by  the  various  executable  files  is  described.  Next  the 
definition  of  control  parameters  and  their  optional  values  is  described.  Finally,  the 
format  and  content  of  the  environment  definition  files  is  described. 

A.1  EXECUTABLE  FILES 

The  M3PE  model  is  executed  by  a  series  of  functions  contained  in  the  files 
listed  in  alphabetical  order  in  Table  5.  All  the  files  in  Table  5  are  executable  by 
MATLAB. 


Table  5.  MATLAB  1 

tile  names  of  functions  executing  the  M3PE  model 

Number 

Filename 

Description 

1 

brough.m 

Computes  a  bottom  roughness  function  that  is  applied 

to  the  mean  bottom  depth.  The  approach  is  based  upon 

Fox  &  Flayes'  definition  as  described  by  Smith, 
Hodgkiss  and  Tappert^"^. 

2 

bw_eff.m 

Computes  the  effective  bandwidth  of  data  in  a  three 

dimensional  matrix,  along  the  last  dimension.  Plots  of 

the  overall  and  conditional  histograms  of  the  features 

are  generated  if  specified. 

3 

get_env.m 

Defines  SSPs,  bottom  depths,  and  environmental 

parameters  vs  range  from  the  options  allowed.  Applies 

desired  bottom  roughness.  Applies  transition  functions 

at  bottom  boundary. 
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4 

init_m3pe.m 

Defines  operating  parameters  that  control  the  execution 

of  the  M3PE  code. 

5 

kfilter.m 

Computes  a  one  sided  taper  function  (sponge)  for 

application  in  k-space. 

6 

mSpe.m 

Main  program  of  the  Matlab  Monterey-Miami  Parabolic 

Equation  model.  Computes  a  transmission  loss  in  range 

and  depth  using  the  Split-Step  Fourier  PE  model. 

7 

m3pe_src.m 

Produces  the  starting  source  field  for  the  M3PE  model. 

It  creates  a  source  field  in  k-f  space  that  accounts  for 

source  motion.  It  also  creates  parameters  that  are  used 

to  implement  effects  of  receiver  motion. 

8 

plot_rkf.m 

Plots  the  wavenumber-frequency  response  at  a 

particular  range. 

9 

raytrace.m 

Traces  a  bundle  of  rays  from  a  source  position,  out  to  a 

maximum  horizontal  range  using  a  Hamiton-Jacobi 

approach,  implemented  with  the  4th  order.  Runge-Kutta 

/  Simpson's  rule  method.  Perfect  reflection  boundaries 

are  included  (if  desired). 

10 

read_env.m 

Reads  the  user  specified  parameters  contained  in  the 

specified  text  file. 

11 

zfilter.m 

Computes  a  one  sided  taper  function  (sponge)  for 

application  in  z-space. 

The  executable  software  also  interfaces  to  two  other  files  used  for  inputting 
acoustical  parameters  and  outputting  the  very  large  three  dimensional  complex 
acoustic  field  computed  when  source  and  receiver  motion  are  greater  than  zero. 
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Table  6.  Input/output  files 


Number 

Filename 

Description 

1 

m3pe_envx.dat 

Contains  parameters  defining  the  range-dependent 

ocean  acoustic  environmental  parameters,  such  as 

water  SSP,  density,  water  absorption  characteristics 

(temperature,  salinity,  pFI),  bottom  sound  speed, 

gradient  of  sound  speed,  density,  and  attenuation. 

Also  defines  the  bottom  depth  versus  range. 

2 

psi.dat 

Temporary  storage  file  of  the  complex  acoustic  field 

vs  range,  depth,  and  frequency. 

The  main  program  is  called  mSpe.m,  and  serves  as  the  origin  for  all  other 
functions.  A  diagram  indicating  the  functional  flow  of  mSpe.m  including  the 
function  dependencies  is  given  in  Figs.  92  and  93.  The  program  is  designed  to 
repeat  its  execution  on  multiple  environments.  The  function  get_env.m  inputs  the 
environmental  data  and  computes  all  the  range-dependent  acoustical 
parameters.  These  are  input  to  the  function  m3pe_src.m  which  computes  the 
starting  field  in  wavenumber  vs  frequency  space.  Taper  functions  applied  to  the 
wavenumber  and  depth  propagators  are  computed  by  the  functions  kfilter.m  and 
zfilter.m,  respectively.  Next,  the  wavenumber  and  depth  propagator  functions  are 
computed  for  the  initial  starting  range.  After  transforming  the  starting  field  from 
wavenumber  to  depth  space,  via  an  inverse  discrete  Fourier  transform  (IDFT), 
the  depth  propagator  is  applied  for  the  first  range  cell.  In  a  loop  on  the  remaining 
range  cells,  the  propagated  field  is  transformed  back  to  k  space  (via  DFT), 
multiplied  by  the  wavenumber  propagator,  and  inverse  transformed  back  to  depth 
space.  Then  the  depth  propagator  is  applied  in  two  stages,  so  that  the  field 
stored  to  disk  can  represent  the  middle  of  a  range  step.  If  source  and/or  receiver 
motion  is  being  modeled,  these  operations  are  repeated  for  all  frequency  bins. 
After  all  range  cells  are  completed,  the  resulting  complex  transmission  loss  is 
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normalized  and  plotted  versus  range,  depth,  and  frequency.  The  magnitudes  of 
the  unnormalized  field  in  three  dimensions  {r,z,f)  are  input  to  the  bw_eff.m 
function  that  computes  aggregate  and  dependent  statistics  on  the  bandwidth  of 
TL  vs  frequency.  Finally,  a  ray  trace  is  computed  for  the  given  SSP  and  bottom 
profile  by  the  function  raytrace.m. 
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Figure  92.  M3PE  Functional  Flow  (part  a) 
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Figure  93.  M3PE  Functional  Flow  (part  b) 
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A.2  CONTROL  PARAMETER  DEFINITIONS 


The  values  of  parameters  that  control  operation  of  the  M3PE  model, 
including  selection  of  various  execution  options,  are  defined  in  the  file 
init_m3pe.m.  This  is  a  MATLAB  executable  file  that  is  included  in  mSpe.m  at  the 
beginning  of  the  file.  A  list  of  all  the  selectable  parameters,  along  with  their 
optional  values  and  nominal  constraints  is  contained  in  Table  7. 


fable  7. 

Control  parameter  definitions 

Num 

Parameter 

Units 

Description 

Normal  Constraints 

1 

f 

Hz 

Transmit  frequency 

5  <  f  <  5000 

2 

zs 

m 

Source  depth 

>  0 

3 

thsd 

degrees 

Source  vertical  angle  of  motion 

-90  <  thsd  <  90 

4 

thrd 

degrees 

Receiver  vertical  angle  of  motion 

-90  <  thrd  <  90 

5 

vsk 

knots 

Source  speed 

-100  <  vsk  <  100 

6 

vrk 

knots 

Receiver  speed 

-100  <  vrk  <  100 

7 

Nf 

integer 

Number  of  frequency  cells 

Nominally  <  100 

8 

nfiles 

integer 

Number  of  environments  to  run 

<  10 

9 

use_file 

Boolean 

integer 

=  1  to  use  a  disk  file  for  temp  storage 

1  or  0 

10 

filename(:,:) 

ascii 

Matrix  of  environment  file  names 

Name  <=  50  chars 

11 

Nr 

integer 

Number  of  range  cells 

32  <=  Nr<=  16384 

(radix  2) 

12 

Nz 

integer 

Number  of  depth  cells 

32  <=  Nz<=  16384 

(radix  2) 

13 

cO 

m/s 

Reference  sound  speed 

1500  m/s 

14 

WARE 

Boolean 

integer 

=  1  to  use  wide  angle  PE  model 

otherwise  use  narrow  angle 

1  or  0 

15 

batten 

Boolean 

integer 

=  1  to  invoke  bottom  attenuation 

1  or  0 
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16 

binterp 

Boolean 

integer 

=  1  to  perform  horizontal  interpolation 

of  bottom  parameters 

1  or  0 

17 

do_cyl_spread 

Boolean 

integer 

=  1  to  include  cylindrical  spreading 

term  in  TL  calculation 

1  or  0 

18 

dorho 

Boolean 

integer 

=  1  to  include  density  in  propagator 

function 

1  or  0 

19 

doalw 

Boolean 

integer 

=  1  to  include  attenuation  in  the  water 

column 

1  or  0 

20 

dorcv 

Boolean 

integer 

=  1  to  accommodate  receiver  motion 

1  or  0 

21 

doraytrace 

Boolean 

integer 

=  1  to  compute  and  plot  a  raytrace 

1  or  0 

22 

dobw 

Boolean 

integer 

=  1  to  compute  bandwidth  statistics 

1  or  0 

23 

dork 

Boolean 

integer 

=  1  to  save  a  k-f  plot  in  the  received 

field 

1  or  0 

24 

maxTL 

dB 

Maximum  TL  for  color  coded  plots 

-500  <  maxTL  <  500 

25 

minTL 

dB 

Minimum  TL  for  color  coded  plots 

-500  <  minTL  <  500 

26 

zplot 

m 

Depth  to  plot  TL  vs  range 

0  <  zplot  <  maxZ 

27 

zbplot 

m 

Depth  to  plot  TL  vs  frequency 

0  <  zbplot  <  maxZ 

28 

Rbplotl 

m 

Min  range  to  plot  TL  vs  frequency 

0  <  Rbplotl  <  maxR 

29 

Rbplot2 

m 

Max  range  to  plot  TL  vs  frequency 

0  <  Rbplot2  <  maxR 

30 

Rkplot 

m 

Range  to  save  k-f  spectral  plot 

0  <  Rkplot  <  maxR 

The  parameters  in  Table  7  are  defined  in  the  file  init_m3pe.m  following  all 
MATLAB  assignment  conventions.  Here  it  is  appropriate  to  discuss  a  few  caveats 
regarding  some  of  these  variables.  First,  all  source  and  receiver  angles  of  motion 
are  positive  down.  Second,  a  positive  source  or  receiver  speed  implies  motion 
towards  increasing  horizontal  coordinates.  Conversely,  a  negative  speed  implies 
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motion  towards  decreasing  horizontal  coordinates.  Regarding  the  mesh  size 
definitions,  it  is  also  important  to  note  that  the  product  of  Nz*Nr*Nf  has  a  limit 
determined  by  the  amount  of  random  axis  memory  (RAM)  installed  in  the 
computer  on  which  the  model  is  being  executed.  The  product  can  be  expanded  if 
the  use_file  option  is  selected,  but  in  either  case,  a  limit  exists  that  is  best 
determined  by  experimentation.  The  nominal  constraints  listed  in  Table  7  were 
selected  based  upon  experience  with  a  computer  containing  512  MB  of  RAM. 

Regarding  the  binary  parameters  listed  in  items  14-23  (in  Table  7),  most  of 
them  are  usually  set  to  1,  when  a  comprehensive  set  of  results  are  desired.  For 
example,  the  WARE  variable  is  only  provided  if  one  desires  to  evaluate  the 
difference  between  wide  and  narrow  angle  approximations  for  academic 
purposes.  The  wide  angle  mode  should  always  be  selected  to  provide  the  best 
accuracy.  Similarly,  the  do_cyl_spread,  batten,  dorho  and  doalw  parameters 
should  only  be  “turned  off”  (i.e.  set  to  zero)  for  special  test  cases.  Obviously 
these  need  to  be  enabled  to  achieve  the  best  accuracy.  The  dorcv  parameter  is 
required  if  the  model  is  to  account  for  receiver  motion.  Next,  the  doraytrace, 
dobw,  and  dork  parameters  simply  control  the  execution  of  post  processing  and 
plotting  functions  and  can  be  selected  according  to  individual  preferences. 

A  cautionary  note  needs  to  be  discussed  regarding  the  binterp  parameter. 
It  should  only  be  enabled  when  the  bottom  acoustical  parameters  are  range 
dependent  but  the  bottom  depth  is  constant.  This  is  because  the  flag  causes  a 
horizontal  interpolation  of  the  bottom  parameters  with  range  that  is  not  yet 
designed  to  handle  variable  bottom  depths.  Finally,  the  parameters  listed  in  rows 
24-30  (of  Table  7)  specify  various  plotting  options  that  are  desirable  for  observing 
the  TL  versus  frequency  response  when  the  source  or  receiver  speeds  are 
greater  than  zero. 


A.3  ENVIRONMENT  DEFINITION  FILES 

This  section  describes  the  format  of  acoustical  parameters  specified  in  the 
ocean  environment  definition  files  referenced  on  row  1  in  Table  6.  These  files 
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allow  specification  of  all  acoustical  parameters  used  by  the  M3PE  model  by 
simply  editing  a  single  ASCII  file.  Although  Table  6  indicates  the  name  of  the  files 
as  having  the  prefix  m3pe_envx,  any  file  name  can  be  assigned  to  the  filename 
variable  listed  in  row  11  in  Table  7.  A  list  of  the  parameters  contained  in  this  file 
is  given  in  Table  8. 


Table  8.  Parameters  contained  in  tl 

he  environmental  definition  files 

Num 

Parameter 

Units 

Normal  Constraints 

1 

Number  of  SSPs 

integer 

>  0 

2 

Range  of  SSP 

m 

0  <  rssp  <  maxR 

3 

SSP  -  depth  vs  sound  speed 

m  -  m/s 

1450<c(z,r)<  1550 

4 

Sound  speed  in  bottom 

m/s 

1450  <c(z,r)<  1550 

5 

Sound  speed  gradient  in  bottom 

1/s 

-20  <  eg  rad  <  20 

6 

Salinity  in  water 

ppm 

30  <  S  <  40 

7 

pH  in  water 

n/a 

1  <  pH  <  20 

8 

Average  temperture  in  water 

°C 

-10<T<50 

9 

Attenuation  in  bottom 

dB/;^ 

0  <  ae  <  10 

10 

Density  profile  in  water  -  depth  vs  rho 

m  -  kg/m^ 

900  <  Pw  <  1500 

11 

Density  in  bottom 

kg/m^ 

1000  <  pw<  3000 

12 

Density  gradient  in  bottom 

kg/m^ 

-20  <  rgrad  <  20 

13 

Bottom  depth  vs  range 

m  -  m 

0  <  Zb  <  10000 

14 

Bottom  roughness  magnitude 

m 

0  <  G  <  100 

15 

Roughness  type 

integer 

1  =  uniform,  2  =  Gaussian,  3  = 

Fox-Hayes,  4  =  sawtooth  right,  5 

=  sawtooth  left 

16 

Correlation  length 

m 

Used  only  by  options  3-5 

Several  of  the  variables  listed  in  Table  7  are  two  dimensional.  For 
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example,  in  rows  3  and  10,  the  sound  speed  and  density  are  defined  as  piece- 
wise  linear  functions  of  depth.  Also,  row  13  refers  to  the  piece-wise  linear 
definition  of  the  bottom  depth  versus  range.  All  other  parameters  are  single 
scalar  variables.  We  also  note  that,  although  the  density  gradient  in  row  12  is 
specified,  it  is  not  yet  used  by  the  model. 

For  the  case  of  a  range-dependent  ocean,  the  parameters  defined  in  rows 
2-12  in  Table  7  are  repeated  for  each  discrete  range  cell.  In  this  case,  the  “range 
of  SSP”  variable  specified  in  row  2  is  assigned  the  range  value  that  corresponds 
to  the  acoustical  parameters  immediately  following  it.  The  “number  of  SSPs” 
variable  defined  in  the  first  row  determines  the  number  of  sets  of  acoustical 
parameters  contained  in  the  file.  After  all  range-dependent  parameter  sets  are 
defined,  the  bottom  depth  profile  and  roughness  parameters  are  listed  at  the  end 
of  the  file. 

A  listing  of  the  contents  of  an  example  file  is  provided  in  Fig.  94.  It  can  be 
observed  that  the  lines  of  ASCII  text  can  be  grouped  into  four  different  types.  The 
first  type  is  preceded  by  a  delimiter  containing  a  single  “%’  character  followed  by 
a  line  of  text.  The  second  type  begins  with  the  “%’  character  and  is  then  followed 
by  a  line  of  dashes  The  third  type  of  line  only  contains  a  pair  of  percent 
characters  (“%%’).  The  final  type  contains  numeric  characters  that  may  or  may 
not  be  followed  by  spaces  and  text  characters. 

The  first  several  lines  in  the  parameter  file,  preceded  by  a  “%’  character, 
can  contain  any  text  entered  by  the  user.  Next  the  parameter  “Number  of  SSPs” 
must  be  preceded  and  followed  by  a  line  of  dashes  as  shown  in  Fig.  94.  After  the 
first  line  of  dashes,  all  the  other  lines  of  text  preceded  by  a  “%’  delimiter  need  to 
contain  the  exact  headings  listed  in  the  figure.  The  user  is  hereby  cautioned  not 
to  modify  these  lines.  For  example  the  text  phrases  “SSP  water  layer,”  “SSP 
bottom  layer,”  and  “Volume  Absorption”  must  all  precede  the  data  corresponding 
to  their  labels.  Also  between  each  type  of  range-dependent  parameters,  e.g. 
“SSP  in  water”,  “Volume  Absorption,”  or  “Bottom  Absorption”  parameters,  the  file 
contains  another  delimiter  represented  by  the  “%%’  characters.  These  are  only 
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inserted  between  the  range-dependent  parameters  that  can  be  repeated 
according  to  the  “Number  of  SSPs”  parameter.  Finally  each  line  containing 
numeric  values  of  some  parameter  can  also  contain  a  series  of  text  symbols  on 
the  same  line.  However  one  or  more  ASCII  spaces  or  tab  characters  must  be 
inserted  between  the  data  and  the  text. 

The  environmental  parameter  file  is  read  by  the  function  read_env.m, 
referenced  in  the  flow  diagram  in  Fig.  92.  This  function  also  has  a  stand-alone 
mode  that  allows  reading  a  file  simply  into  the  MATLAB  workspace.  Then  the 
variables  can  be  individually  evaluated  to  ensure  consistency  with  expectations. 
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%  Fie  m3pe_env1  .dat 
%  This  He  contans  parameters  that  describe 
%  the  acoustic  envionment  of  a  modeled  ocean 
% - 

%  Number  of  SSPs 
1 

% - 

%  SSP  vrater  layer  (to  max  bot  dep  spac'd  below) 
0  range  associated  with 

%  depth  (m)  sound  speed  (m/s) 

0  1500 

200  1480 

3000  1540 

%% 

%  SSP  bottom  layer  parameters 
1600  Sound  speed 

.2  Gradient  (del  SOS  /  del  depth) 

%% 

%  Vohime  Absorption  parameters 
35  Sairiity  (ppm) 

8  pH 

15.6  Tempertiae  (deg  C) 

%% 

%  Bottom  Absorption  parameter 
0.1  dB/bmbda 

%% 

%  DP  of  vrater  (Dertsity  vs  depth  of  vrater)  layer 
%  depth  (m)  dens%  (kg/m*3) 

0  1024 

3000  1024 

%% 

%  DP  of  bottom  (Dens%  params  of  bottom)  layer 
1700  Dens% 

0  Gradient  (del  dens  /  del  depth) 

% - 


%  Bottom  depth  vs  range  (m) 

%  depth  (m)  range  (m) 

3000  0 

3000  10000 

2000  20000 

1000  30000 

500  35000 

500  40000 

% - 

%  Bottom  toughness  parameter 
1  roughness  pk-pk  magnitude  (m) 

3  rtype,  1  =  uniform,  2  =  Gaussian,  3  =  Fox- 

Hayes,  4  =  sawtooth  right  5  =  sawtooth  left 
100  Lst  =  length  of  sawtooth  (m)  if  appkcable 
% - 


Range  dependent. 
Can  be  repeated 
according  to  the 
Number  of  SSPs 
definition  at  the  top. 


Figure  94.  Example  ocean  environment  definition  file 
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APPENDIX  B.  RANGE  DEPENDENT  RAY  TRACE  FUNCTION 


This  appendix  presents  the  derivation  of  the  method  used  to  produce  the 
ray  traces  included  throughout  this  thesis.  The  first  section  derives  the 
mathematics  associated  with  ray  theory,  and  the  second  section  describes  the 
numerical  technique  which  solves  the  equations.  A  third  section  provides  a 
description  of  the  MATLAB  code  that  implements  the  functions. 

B.1  RAY  THEORY 

The  first  part  of  our  discussion  below  closely  follows  parts  of  Refs.  41  and 
42.  We  begin  by  restating  the  homogenous  Helmholtz  equation  as 

VV  +  -^/’  =  0  ,  (B1) 

c  (x) 

where  x  is  the  cartesian  space  coordinate.  As  usual,  we  can  again  assume  that 
the  pressure  is  a  spatially  variant  time  harmonic  function  defined  by 

p{x)=  .  (B2) 

Then  the  spatial  gradient  is 

Vp  =  e“^^(VA  +  i(i)AVr)  ,  (B3) 

and  the  Laplacian  is 

Vp^  =  e“"^{ia)VA  •  Vr  -  +  V^A  +  ioNA  •  Vr  +  ia)AV\)  .  (B4) 

Substituting  (B4)  into  (B1),  and  dropping  the  time-dependent  term,  results  in 

iQNA-VT-Q/A\VT\  +V^A  +  i(i>VA-VT  +  i(i>AV^T  +  —  A  =  0  .  (B5) 

Equation  (B5)  is  true  when  both  the  real  and  imaginary  parts  are  equal  to 
zero.  Thus, 
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(B6) 


V^A-co^AIVtI^  +^A  =  0 
c 


and 


(yVy4  •  Vr  +  (yVy4  •  V  r  +  =  0  .  (B7) 

Next  we  take  the  limit  of  Eqs.  (B6)  and  (B7)  as  ^  oo ,  which  results  in 

|Vr|-=i  (B8) 


and 


2V^-Vr  +  ^V\  =  0  .  (B9) 

Equation  (B8)  is  commonly  called  the  eikonal  equation,  while  (B9)  is  called  the 
transport  equation.  The  eikonal  equation  determines  the  trajectory  of  vectors  that 
are  perpendicular  to  the  wavefront,  or  position  of  constant  phase.  The  transport 
equation  determines  the  amplitude  along  the  wavefront  trajectory.  To  derive  the 
ray  tracing  formulae,  we  focus  on  the  eikonal  equation. 

Equation  (B8)  can  be  rewritten  in  cylindrical  coordinates  as 


dr 

\Srj 


r  1 


+ 


ydzj 


2  r 
+ 


1  dr 
dcj) 


(B10) 


Here  we  note  that  the  speed  of  sound  is  defined  to  be  generally  variable  with  all 
spatial  coordinates.  If  we  neglect  azimuthal  variations,  this  becomes 


ydr  j 


ydz  J 


(B11) 


Eqs.  (B10)  and  (B11)  now  have  the  form  of  a  Hamilton-Jacobi  equation 
used  to  describe  the  classical  motion  of  a  particle.  The  analogy  between  sound 
propagation  and  particle  motion  has  been  known  for  many  years  (see  [43]  for 
example).  The  following  discussion  shows  how  methods  from  classical 
mechanics  provide  a  convenient  framework  for  solving  the  acoustic  wavefront 
trajectory  problem. 
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Let  us  briefly  digress  to  highlight  some  of  the  aspects  of  this  analogous 
theory.  Our  discussion  closely  follows  the  contents  of  reference  [44],  Hamilton’s 
principle  states  that  the  actual  path  taken  by  a  particle  between  fixed  end  points 
will  be  the  path  that  minimizes  a  quantity  called  the  “action”.  The  action  is  defined 
as  the  integral  of  the  difference  between  the  kinetic  energy,  T  (a  function  of 
velocity),  and  the  potential  energy,  V  (a  function  of  position),  defined  according  to 

tl 

5  =  j(r(4-F(x))Jr  .  (B12) 


Here  the  variable  x  defines  the  position,  and  ti  and  t2  are  the  time  end  points  for 
which  the  particle  motion  is  being  evaluated.  The  integrand  is  also  called  the 
Lagrangian,  defined  as  l{x,:^=  t{j^-V{x) .  The  variable  &  refers  to  the  time 
rate  of  change  ofx.  In  Hamiltonian  mechanics,  the  momentum  associated  with  a 
generalized  coordinate,  qj,  is  a  function  of  the  Lagrangian  according  to 


Pi 


(B13) 


The  momentum  p  .,  the  generalized  coordinate  variable  qj,  and  the  Lagrangian, 
L{q,<^,  together  define  the  Hamiltonian  function,  H{p,q)  ,  as 


!  =  1 


(B14) 


where  D  is  the  dimensionality  of  the  generalized  coordinates,  i.e.  the  number  of 
degrees  of  freedom  in  the  system. 

The  Hamiltonian  is  analogous  to  the  total  energy  of  the  system"^^.  It  is 
related  to  the  generalized  coordinate  and  the  momentum  by  the  following  first 
order  equations: 


Pj  = 


dH 


(B15) 


117 


dH 

dpj 


(B16) 


Equations  (B15)  and  (B16)  are  referred  to  as  Hamilton’s  equations  of  motion. 
They  represent  a  pair  of  first-order  coupled  differential  equations  that  can  be 
more  easily  solved  than  a  higher  order  set  of  equations. 


Given  these  tools,  and  recognizing  r  as  the  action-like  variable  and  r  as 
the  time-like  variable  in  Eq.  (B11),  we  can  begin  to  solve  (B11)  by  defining  the 
momentum  term,/>,  as 


P 


dz 


(B17) 


Then,  Eq.  (B11)  can  be  rewritten  as 


dr 

dr 


(B18) 


is 


Next,  we  can  define  the  term 


Qt 

—  as  the  negative  of  the  Hamiltonian.  That 

dr 


H{p,z,r)  =  -^  =  - 
dr 


ir,zY 


(B19) 


Then,  by  using  (B15)  and  (B16),  it  follows  that 

dp  dH 


dr  dz 


(B20) 


and 


dz  dH 
2^ —  — 

dr  dp 


(B21) 


Substituting  (B19)  into  (B20)  and  taking  the  derivative  results  in 


118 


dp 

dr 


-1 


dc{r,  z)  _ 


-1 


dc{r,  z) 


ir,zy 


ir,zf 


dz 


c{r,zy{^-c{r,zY 


dz 


Performing  this  same  substitution  into  Eq.  (B21)  also  results  in 

p  _  c{r,z)p 


dz  _ 
dr  / 


1 


:(r,z)' 


2  (l-c(r,z)V7 


(B22) 


(B23) 


Now  recall  that  —  is  nothing  more  than  the  slope  of  the  ray  perpendicular 

dr 

to  the  direction  of  the  propagating  wavefront  in  range  and  depth  space.  The 
relationship  of  this  derivative  to  the  grazing  angle  of  propagation  is 

tan(6»)  =  —  .  (B24) 

dr 


Thus, 


—  =  tan(^)  = 
dr 


sin(^)  _  c{r,z)p 


(B25) 


Equating  the  numerators  on  both  sides  of  (B25)  results  in  the  following  simple 
definition  for  the  momentum: 


P  = 


1 

c(r,z) 


sin(6’)  . 


(B26) 


In  (B22),  (B23),  and  (B26),  we  now  have  a  system  of  equations  that  can 
be  numerically  integrated  to  solve  for  the  trajectory  in  range  and  depth  of  a 
perpendicular  ray  to  an  acoustic  phase  front.  Prior  to  discussing  numerical 
methods,  we  also  need  an  expression  for  the  travel-time  of  a  ray  through  space. 
This  can  be  attained  by  using  a  method  from  classical  mechanics  as  well.  Again 
we  closely  follow  reference  [42].  In  this  case,  Fermat’s  Theorem  tells  us  that  the 
travel  time  of  a  particle  is  related  to  the  Lagrangian  by 
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rl 


T  =  ^L{r,z)dr  . 


From  (B14)  the  Lagrangian  is  related  to  the  Hamiltonian  by 
L{z,r)  =  ^-H{p,z,r)  . 

Also  recall  the  definition  of  the  Hamiltonian  as 


H{p,z,r)  =  - 


1 


ir,zf 


J 


-^{^-c{r,zY  p^\ 
ir,z) 


Substituting  (B23)  and  (B29)  into  (B28)  results  in 

c{r,zfp^  ,  i^-c{r,zY  p^) 


L{z,  r)  = 


-  +  - 


ir,  z)(l  -  c(r,  zY  p^y  c{r,  z)(l  -  c(r,  z)'  p "  Y 

This  can  be  simplified  to 

L{z,r)=  ^  ^ 


1 


(B27) 


(B28) 


(B29) 


(B30) 


ir,z){l-c{r,zYp^y  .  /  .  . 

c(r,z)  l-c(r,z) 


■  (B31) 


ir,zY 


Thus  the  integral  of  (B31)  along  the  ray  path  provides  the  travel  time  of  the  ray 
between  the  endpoints  of  integration.  That  is 


r2 


(r,z)cos(^) 


dr  . 


(B32) 


B.2  NUMERICAL  SOLUTION 

We  now  describe  a  procedure  for  generating  the  range  and  depth 
coordinates  of  an  acoustical  ray  originating  from  a  source  position  with  a 
particular  vertical  launch  angle.  Our  problem  requires  integration  of  Eqs.  (B22) 
and  (B23)  over  a  horizontal  range  interval.  This  is  performed  using  the  Runge- 
Kutta  Method,  as  described  by  Smith"^^"^^  for  this  application.  Our  discussion  also 
relies  heavily  upon  reference  [46]. 
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Equations  (B22)  and  (B23)  follow  the  general  form  of  a  coupled  first  order 
ordinary  differential  equation  that  can  be  re-written  as 

^  =  f{x,y{x))  .  (B33) 

ax 

This  implies  that  the  differential  of  a  function,  y(x),  is  dependent  upon  the  original 
function,  and  an  independent  variable,  x.  Problems  of  this  sort  can  be  solved  by 
the  recursion"^"^, 

Xn  +/? 

T«+i  =T«  +  >  (B34) 

where  h  is  the  length  of  the  sub-interval  along  the  independent  variable  x.  In  the 
4*^^  order  Runge-Kutta  approach,  the  solution  is  obtained  numerically  by 


7  (^1  +  2^2  +  2^3  +  ^4 )  (B35) 

6 

where 

ki=h-f{x„,yj  ,  (B36a) 

ki^h- f(x„+^h,y„+^k,^  ,  (B36b) 

k,=h-  f(x„+^h,y„+^k2]  ,  (B36c) 

k,=h-f{x„+h,y„+k,)  .  (B36d) 


It  is  also  interesting  to  note  that  if  the  function/)  is  independent  of  j(x),  then  k2  = 
ks  and  (B36)  reduces  to  Simpson’s  rule,  (where  h  is  substituted  for  h/2),  i.e., 

f{x„)  +  ^f[x„+^h^  + f{x„+h)  ■ 

We  are  now  ready  to  show  how  the  above  methods  can  be  applied  to  our 
problem.  For  completeness,  we  restate  Eqs.  (B22)  and  (B23)  below: 
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dp  _  - 1 


dr 


c{r,zY{^-c{r,zY 


dc{r,z) 

dz 


dz 


c{r,z)p 


dr 


^-c{r,zY  p^y 


(B22) 


(B23) 


In  our  case,  the  range  variable  r  plays  the  role  of  the  independent  variable  x  in 
(B36).  Also,  the  range  increment  Ar  plays  the  role  of  the  interval  length  h.  The 
recursion  relation  for  depth  z  can  then  be  found  from"^^ 


^«+l  -^n  +7(^1  +2^2  +2^:3  +k^) 
6 


(B37) 


K=/kr^{z„,Pn,r)  , 
dr 


(B38a) 


^2 


r  1  ,  1 

^n+i:k^,p„+-m„r  +  —  , 

\  Z  Z  Z  J 

f  \  ,  1  Ar'] 

+:zk2,Pn  +i:m^^r  +  — 

\  Z  Z  Z  ) 


(B38b) 

(B38c) 


k. 


^k^,p^  +m„r 
dr 


+  Ar) 


(B38d) 


A  similar  set  of  steps  must  be  simultaneously  executed  for  the  momentum,  p 
according  to 


Pn^i  =  ;?„  +  7(^1  +  2^2  +2m^+m,)  , 
6 


(B39) 


My 


^r^{Pn,z„,r) 

dr 


(B40a) 


A  dp 

m-,  =  Ar  — 
"  dr 


=  Ar 


dp 

dr 


1 


1 


Ar 


o„  H — m,  ,z„  H — k,,r  - 

2  ‘  "  2  ‘  2 


1  1  ,  Ar 

Pn  +^k„r  +  — 


(B40b) 

(B40c) 
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^r^{p,+m^,z„+k^,r  +  Ar)  . 
dr 


(B40d) 


= 


Our  final  step  is  to  present  the  approach  for  computing  the  travel  times  by 
integrating  the  Lagrangian  as  presented  in  (B32).  Smith"^^  recommends  using 
Simpson’s  rule, 

r  =  ^[z(r)  + 4Z(r  +  Ar)  + 2Z(r  +  2Ar)+ ...  +  4Z(r  +  (n  -  l)Ar)  + Z(r  +  nAr)]  .  (B41) 

Here  n  is  the  number  of  range  cells  for  which  the  Lagrangian  (specified  by  (B31)) 
is  evaluated.  Thus,  the  Lagrangian,  A(r)=  [c(z,r)cos(6’(z,r))]  ' ,  is  evaluated  at 
each  position  along  the  ray  path  and  plugged  into  (B41)  to  compute  the  travel 
time  up  to  a  particular  position  along  the  ray.  The  propagation  angle,  0,  can  be 
found  from  (B26),  i.e. 

6{z,r)  =  ?,var^\c{z,r)p{z,r^  ’  .  (B42) 


B.3  MATLAB  IMPLEMENTATION 

A  MATLAB  function  has  been  written  to  implement  the  method  described 
in  the  previous  sub-section.  The  function  is  called  raytrace,  and  is  executed  by 
calling  the  function  inside  a  main  program.  The  inputs  to  the  function  are  listed 
and  described  in  Table  9.  These  inputs  are  specified  in  the  argument  list  of  the 
raytrace. m  function. 

It  is  noted  that  the  sound  speed  and  bottom  depth  can  be  range- 
dependent.  However,  due  to  the  simplicity  of  the  reflection  method,  it  is  not 
recommended  that  bottom  depths  with  applied  roughness  be  input  to  this  routine. 
In  the  M3PE  implementation,  a  mean  bottom  depth  (without  roughness)  was 
used  to  create  all  the  ray  trace  plots  illustrated  in  the  various  sections  of  this 
thesis. 
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Table  9.  Input  Variables  to  the  Raytrace  Function 


Number 

Variable 

Description 

1 

zs 

Source  depth 

2 

z 

Depth  vector 

3 

Nz 

Number  of  depth  cells 

4 

c 

Speed  of  sound  matrix  (range  vs  depth) 

5 

maxR 

Maximum  range  (m) 

6 

zb 

Bottom  depth  vs  range 

7 

itheta 

Center  angle  in  degrees  of  rays  at  the  source 

8 

dtheta 

Width  of  ray  bundle  in  degrees 

9 

delth 

Angle  separation  of  rays  at  the  source 

10 

nobounce 

=  1  to  suppress  boundary  interactions 

11 

pssp 

=  1  to  plot  the  ssp,  0  otherwise 

12 

istatus 

=  1  to  plot  status  messages 

The  outputs  of  the  function  consist  of  color  plots  of  the  ray  trace  and  the 
travel  time  versus  range.  Each  ray  is  given  a  separate  color  or  line  characteristic 
(solid  versus  dashed)  to  facilitate  following  the  ray  through  range  and  depth 
space.  A  maximum  of  12  separate  line  types  are  available.  When  greater  than  12 
rays  are  plotted,  the  color  or  line  type  characteristics  are  reused  in  modulo 
fashion.  Optionally,  the  sound  speed  profile  can  also  be  plotted.  Finally,  the  travel 
times  versus  launch  angle  and  depth  at  the  final  range  are  printed  to  the 
workspace  window  in  a  formatted  list. 

The  processing  sequence  within  the  raytrace. m  routine  is  illustrated  by  the 
flow  charts  in  Figs.  97  and  98.  It  can  be  observed  that  the  numerical  integration  is 
performed  inside  two  nested  loops,  over  the  range  grid  and  for  each  individual 
ray. 

When  a  boundary  interaction  is  recognized,  a  series  of  linear  interpolation 
operations  are  executed  depending  upon  whether  the  boundary  is  the  bottom  or 
the  surface.  Figure  95  is  provided  to  help  illustrate  our  notation  for  the  bottom 
interaction.  In  this  case  the  incident  ray  impinges  on  the  bottom  at  an 
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angle  of  6o  relative  to  the  horizontal  axis  (positive  angles  down).  The  last  sample 
of  the  ray  before  the  bottom  is  exceeded  is  located  at  [rn-i,  where  “n”  refers 
to  the  range  index. 


Figure  95.  Bottom  reflection  geometry 


The  local  angle  of  the  bottom  slope  is 


6*^  =  tan  ' 


f  b  b  \ 
Z  —  Z  , 

n  n-\ 

Ar 


(B43) 


where  z*  is  the  bottom  depth  at  the  discrete  range  increment  “n”,  and  Ar  is  the 
range  resolution.  Then  the  reflected  momentum  is 


p  =  sin 

V 


0^ 

c{z,r)) 


(B44) 


Next,  the  depth  of  the  ray  after  the  reflection  is  computed  for  the  discrete  range 
associated  with  index  “n”  by  the  following  steps.  The  angle  between  the  incident 
ray  and  the  bottom  is 
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(B45) 


The  angle  between  the  reflected  ray  and  the  horizontal  axis  is 

6' =26,-6,  .  (B46) 


The  equation  of  the  line  defining  the  incident  ray  is 


T 


Ar 


x  +  z, 


n-\  ’ 


(B47) 


and  the  equation  of  the  local  bottom  is 


T  = 


Ar 


X  +  z 


b 

n-\ 


(B48) 


Equating  vertical  variables  provides  a  solution  to  the  horizontal  component  of  the 
incident  ray, 


Ax  = 


(B49) 


Here  the  origin,  z*,  is  the  point  at  which  the  ray  intersects  the  bottom.  The 
horizontal  component  of  the  reflected  ray  is 

Ax'=Ar-Ax  .  (B50) 

The  vertical  component  of  the  reflected  ray  is 

Ay=  Ax'tan(^’)  .  (B52) 

Finally,  the  new  depth  of  the  ray  relative  to  the  surface  is 

K=4-^y'  ■  (B53) 

After  the  reflection,  the  ray  is  propagated  from  the  point  [r^,z„]  with  the  new 
momentum  defined  by  Eq.  (B44), 

If  the  ray  depth  during  the  propagation  loop  becomes  less  than  zero  (at 
the  surface)  a  similar  approach  is  followed.  An  illustration  of  the  geometry  is 
provided  in  Fig.  96. 
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Figure  96.  Surface  reflection  geometry 


The  ray  depth  after  reflection  is  computed  from  the  following  steps.  If  0^  is  the 
angle  of  the  incident  ray  (positive  down),  then 

e,  =  -e,  (B54) 

is  the  angle  between  the  incident  ray  and  the  horizontal  axis.  Then  the  horizontal 
component  of  the  ray  is 


Ax  = 


tan(6*j) 


The  horizontal  component  of  the  reflected  ray  is 

Ax'=Ar-Ax  . 


Finally  the  vertical  component  is 

z„  =  Ax'  tan(6'i )  . 

From  the  point  [r„,z^],  the  ray  is  propagated  with  the  new  momentum 


(B55) 


(B56) 


(B57) 


p  =  sin 

V 


c{z,r)) 


(B58) 
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Figure  97.  Flow  Diagram  of  Raytrace  Function,  part  (a) 
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Figure  98. 
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Flow  Diagram  of  Raytrace  Function,  part  (b) 
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B.4  BENCHMARK  EXAMPLES 


To  illustrate  its  capabilities,  the  raytrace  program  was  executed  with 
Munk’s  canonical  SSP  down  to  a  depth  of  4500  m.  A  flat  bottom  was  assumed 
out  to  a  range  of  200,000  m.  The  SSP  is  explicitly  defined  as  follows: 


c 


(^) 


=  Cn 


(B43) 


2(z  —  z  .  ) 

where  tj  = - and  =  1000  m ,  5  =  1000  m,  ^  =  0.0057  ,  and 

B 

Cq  =  1490  m/s  .  A  plot  of  this  SSP  is  illustrated  in  Fig.  99. 


Speed  of  Sound  (m/s) 


Figure  99.  Munk’s  canonical  profile 


A  mesh  size  of  Az  =  400  and  N,-  =  8000  was  selected  for  our  plot.  With  the 
source  depth,  zs  (from  Table  9)  equal  to  1000  m,  the  rays  propagating  from 
launch  angles  between  [-30°,30°]  at  0.5°  increments  are  plotted  in  Fig.  100. 
Figure  101  illustrates  the  trace  from  a  source  at  500  m.  Reference  41  can  be 
consulted  to  find  examples  using  the  Munk  profile  and  a  flat  bottom  for 
comparison. 
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Ray  Trace  Plot,  60  rays,  start  angles  [-15.0,14.5]  degs,  zs  =  1000  m 
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Figure  100.  Example  ray  trace,  Zg  =  1000  m 


Figure  101.  Example  ray  trace,  Zg  =  500  m 
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APPENDIX  C.  SCALAR  FORM  OF  THE  KE  OPERATOR 


This  appendix  derives  a  scalar  version  of  the  KE  operator  defined  in  Eq. 
(20)  that  is  equivalent  to  the  operator  form  given  by  Eq.  (16).  The  result  is  used 
as  a  multiplier  in  the  wavenumber  domain  within  the  Split-Step  Fourier  marching 
algorithm  defined  by  Eq.  (21).  We  begin  by  restating  Eq.  (16)  as 


Top  =\-  +  M  =  1 


(16) 


Equation  (16)  can  be  generalized  as 


fop{x)  =  \ 


(Cl) 


Then  fop{x)  can  be  expanded  into  a  Maclaurin  series  as 


^OP  (•^)  “  ^ 


A. 


0  n\ 


(C2) 


where  the  coefficients  An  are  the  successive  derivatives  of  T^p(x)  evaluated  at 
zero  according  to 

l  =  •  (C3) 

Also 


X 


1  5 


dz 


(C4) 


Similarly,  the  operator  Top  can  be  rewritten  as 


'-OP 


V&y  „=o 


Ail 

n\  dz" 


(C5) 


where  the  exponential  factor  has  been  incorporated  directly  into  the  partial 
derivative. 
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Recall  from  Eq.  (15)  that  Top  is  applied  to  the  PE  field  function  \)/(z).  It  then 
follows  that 


*-OP 


ydzj 


A..  5” 


!-(-)= • 

=n  n\  oz 


(C6) 


«=0 


Then  \|/(z)  can  be  written  in  terms  of  its  wavenumber  transform  as 


y/(z)  =  ^^(k^)e‘^^^dk^ 


(C7) 


Substituting  Eq.  (C7)  into  (C6)  results  in 


'-OP 


\dzj 


(C8) 


In  (C8),  it  is  possible  to  interchange  the  partial  derivative  and  integral  since  only 
the  exponential  term  in  the  integral  is  a  function  of  z.  The  result  is 


'-OP 


'ydzj 


n=o  r^z 


ez" 


After  taking  the  derivatives,  the  result  is 


'^OP 


\dzj 


¥{z)  =  Z^I  A:‘{kX^ky  e’^^^dk^ 

n=0  n\ 


(C9) 


(CIO) 


Then  all  the  functions  of  the  variable  n,  including  the  summation,  can  be  brought 
inside  the  integral  as 


'OP 


ydzj 


^dk. 


(C11) 


Now  the  terms  in  the  parenthesis  in  the  right  side  of  Eq.  (C11)  have  the  same 
form  as  Eq.  (C2).  Therefore  we  can  rewrite  (Cl  1 )  as 


'OP 


ydzj 


y/{z)  =  J  A>{k^)fop{k^)e'‘^^dk^ 


(C12) 


Since  Eq.  (C2)  is  equivalent  to  Eq.  (Cl),  it  follows  that 
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(C13) 


Thus,  Eq.  (C12)  tells  us  that  applying  the  operator  Top  to  the  field  function  \|/(z)  is 
equivalent  to  multiplying  the  scalar  function  fQp{k^)  by  the  transformed  field 
and  transforming  the  result  back  into  the  z  domain. 
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